Options and VaR with R
Dr. Martín Lozano.
Current 18 dic. 2022, 15:05:48.
Abstract
This material relies on John C. Hull (Hull 2015). Some mathematical background is skipped to emphasize the data analysis, model logic, discussion, graphical approach and R coding (R Core Team 2023). As in the philosophy of Donald Knuth (Knuth 1984), the objective of this document is to explain to human beings what we want a computer to do as literate programming. This is a work in progress and it is under revision.
1 Payoff functions.
The payoff is the cash realized by the holder of an option or other asset at the end of its life.
1.1 European options.
It is often useful to characterize a European option in terms of its payoff to the purchaser of the option. The initial cost of the option is then not included in the calculation. The payoff functions for European options are the following.
- Payoff from a long position in a European call option: \(\mathrm{max}(S_T-K, 0)\). Or equivalently: \(-\mathrm{min}(K-S_T, 0)\).
- Payoff to the holder of a short position in the European call option: \(\mathrm{min}(K-S_T, 0)\). Or equivalently: \(-\mathrm{max}(S_T-K, 0)\).
- Payoff to the holder of a long position in a European put option: \(\mathrm{max}(K-S_T, 0)\). Or equivalently: \(-\mathrm{min}(S_T-K, 0)\).
- Payoff from a short position in a European put option: \(\mathrm{min}(S_T-K, 0)\). Or equivalently: \(-\mathrm{max}(K-S_T, 0)\).
Options are zero sum games because the payoffs to the holder or a long position and to the holder of a short position offset each other.
In the case of a call option:
\(\mathrm{max}(S_T-K, 0) + \mathrm{min}(K-S_T, 0)\),
\(\mathrm{max}(S_T-K, 0) - \mathrm{max}(S_T-K, 0) = 0\).
In the case of a put option:
\(\mathrm{max}(K-S_T, 0) + \mathrm{min}(S_T-K, 0)\),
\(\mathrm{max}(K-S_T, 0) - \mathrm{max}(K-S_T, 0) = 0\).
ST.seq <- seq(from = 0, to = 150, length.out = 150)
K = 75
cl <- pmax(ST.seq - K, 0) # long call payoff.
pl <- pmax(K - ST.seq, 0) # long put payoff.
cs <- pmin(K - ST.seq, 0) # short call payoff.
ps <- pmin(ST.seq - K, 0) # short put payoff.Graphically (assuming \(K=75\)):
par(mfrow = c(2, 2), mai = c(0.7, 0.4, 0.4, 0.4))
par(pty = "s")
plot(ST.seq, cl, type = "l", ylim = c(0, 150), lwd = 4,
ylab = "Payoff", xlab = "")
legend("topleft", legend = c("Long call"),
bg = "white", bty = "n")
par(pty = "s")
plot(ST.seq, pl, type = "l", ylim = c(0, 150), lwd = 4,
ylab = "", xlab = "")
legend("topleft", legend = c("Long put"),
bg = "white", bty = "n")
par(pty = "s")
plot(ST.seq, cs, type = "l", ylim = c(-150, 0), lwd = 4,
ylab = "Payoff", xlab = expression(paste(S[T])))
legend("bottomleft", legend = c("Short call"),
bg = "white", bty = "n")
par(pty = "s")
plot(ST.seq, ps, type = "l", ylim = c(-150, 0), lwd = 4,
ylab = "", xlab = expression(paste(S[T])))
legend("bottomleft", legend = c("Short put"),
bg = "white", bty = "n")
Figure 1.1: Payoffs of European option positions.
1.2 Stocks.
The payoff functions for stocks are the following.
Long position in a stock: \(S_T-S_0\).
Short position in a stock: \(S_0-S_T\).
Graphically:
S0 = 75
par(pty = "s")
par(mfrow = c(1, 2), oma = c(0, 0, 2, 0))
par(pty = "s")
# Long position
plot(ST.seq, ST.seq-S0, type = "l", ylab = "Payoff", main = "Long stock",
xlab = expression(paste(S[T])), lwd = 4, col = "blue",
ylim = c(-75, 75))
abline(h = 0, lty = 2)
abline(v = 0, lty = 2)
par(pty = "s")
# Short position
plot(ST.seq, S0-ST.seq, type = "l", ylab = "", main = "Short stock",
xlab = expression(paste(S[T])), lwd = 4, col = "red",
ylim = c(-75, 75))
abline(h = 0, lty = 2)
abline(v = 0, lty = 2)
Figure 1.2: Long and short stock payoffs.
1.3 Bonds.
The payoff functions for zero coupon bonds are the following.
Long position in a bond: \(K\).
Short position in a bond: \(-K\).
Graphically (assuming \(K=75\)):
K = 75
par(pty = "s")
par(mfrow = c(1, 2), oma = c(0, 0, 2, 0))
par(pty = "s")
# Long position
plot(ST.seq, rep(K, 150), type = "l", ylab = "Payoff", main = "Long bond",
xlab = expression(paste(S[T])), lwd = 4, col = "blue",
ylim = c(0, 150))
abline(h = 0, lty = 2)
abline(v = 0, lty = 2)
par(pty = "s")
# Short position
plot(ST.seq, rep(-K, 150), type = "l", ylab = "", main = "Short bond",
xlab = expression(paste(S[T])), lwd = 4, col = "red",
ylim = c(-150, 0))
abline(h = 0, lty = 2)
abline(v = 0, lty = 2)
Figure 1.3: Long and short zero coupon bond payoffs.
1.4 Options, stocks and bonds.
We look at what can be achieved when an option is traded in conjunction with other assets. In particular, we examine the properties of portfolios consisting of (a) an option and a zero-coupon bond, (b) an option and the asset underlying the option, and (c) two or more options on the same asset. This section is based in Hull (2015) chapter 12, you can also see Brealey et al. (2020).
S0 = 65
# Here I assume the stock price is currently 65.
Sl <- ST.seq - S0
# I assume the option price is 5.
cs.pr <- pmin(K - ST.seq, 0) + 5
# A short put profit is replicated by taking a short call and a long stock.
ans <- cs.pr + Sl
par(pty = "s")
plot(ST.seq, ans, type = "l", ylim = c(-75, 75), lwd = 2, col = "red",
xlab = expression(paste(S[T])),
ylab = "Profit")
lines(ST.seq, Sl, lty = 2, lwd = 2, col = "orange")
lines(ST.seq, cs.pr, lty = 2, lwd = 2, col = "blue")
abline(v = 0)
abline(h = 0)
abline(v = 65, lty = 2, col = "grey")
abline(v = 75, lty = 2, col = "grey")
abline(v = 80, lty = 2, col = "grey")
legend("topleft", legend = c("Short call", "Long stock", "Short put"),
col = c("blue", "orange", "red"), lty = c(2, 2, 1),
lwd = 2, bg = NULL)
Figure 1.4: Hull: Figure 12.1a, profit diagram. A synthetic short put.
S0 = 65
c_price = 5
# Here I assume the stock price is currently 65.
Sl <- ST.seq - S0
# I assume the option price is 5.
cs.pr <- pmin(K - ST.seq, 0) + c_price
# A short put profit is replicated by taking a short call and a long stock.
ans <- cs.pr + Sl
par(pty = "s")
plot(ST.seq, ans, type = "l", ylim = c(-5, 15), xlim = c(60, 80),
lwd = 4, col = "red",
xlab = expression(paste(S[T])),
ylab = "Profit")
lines(ST.seq, Sl, lwd = 4, col = "orange")
lines(ST.seq, cs.pr, lwd = 4, col = "blue")
abline(v = 0)
abline(h = 0)
abline(v = 65, lty = 2, col = "grey")
abline(v = 75, lty = 2, col = "grey")
abline(v = 80, lty = 2, col = "grey")
points(S0 - c_price, 0, cex = 3, col = "red", lwd = 3)
points(K, K - S0 + c_price, cex = 3, col = "red", lwd = 3, pch = 0)
points(K + c_price, 0, cex = 3, col = "blue", lwd = 3)
points(K, c_price, cex = 3, col = "blue", lwd = 3, pch = 0)
points(S0, 0, cex = 3, col = "orange", lwd = 3)
points(K, K - S0, cex = 3, col = "orange", lwd = 3, pch = 0)
points(K, 0, cex = 3, lwd = 3)
Figure 1.5: Hull: Figure 12.1a, profit diagram (zoom). A synthetic short put.
We can identify at least three isosceles right triangles, they have two sides of equal measure. The short call profit is zero at \(S_T = 80\), in the code above: K + c_price. The stock profit is zero at \(S_T=65\), in the coded above: S0. Finally, the short put profit is zero at \(S_T=60\), in the code above: S0 - c_price.
Now see how it looks if we plot a payoff diagram instead.
S0 = K
Sl <- ST.seq - S0
cs.pr <- pmin(K - ST.seq, 0)
# A short put profit is replicated by taking a short call and a long stock.
ans <- cs.pr + Sl
par(pty = "s")
plot(ST.seq, ans, type = "l", ylim = c(-75, 75), lwd = 5, col = "red",
xlab = expression(paste(S[T])),
ylab = "Payoff")
lines(ST.seq, Sl, lty = 2, lwd = 2, col = "orange")
lines(ST.seq, cs.pr, lty = 2, lwd = 2, col = "blue")
abline(v = 0, lty = 2)
abline(h = 0, lty = 2)
legend("topleft", legend = c("Short call", "Long stock", "Short put"),
col = c("blue", "orange", "red"), lty = c(2, 2, 1),
lwd = 2, bg = "white")
Figure 1.6: Hull: Figure 12.1a, payoff diagram.
library(knitr)
library(kableExtra)
library(plotly)
library(ggplot2)
library(tidyr)
library(dplyr)
library(vembedr)2 Options properties.
Financial options are one of the most flexible, popular and versatile financial instruments. They are relevant by themselves and important as building blocks to create other financial products and financial strategies. In order to fully understand options we need to understand their properties first. This is, what are the determinant of financial options, how these determinants change the option prices, what are the limits of financial options, what happens when these limits are violated, methods for determining the theoretical option prices and the relationship between call and put options. This section is based in Hull (2015) chapter 11.
2.1 Determinants.
The main determinants of option prices are those parameters of the Black-Scholes formula:
- The price of the underlying stock at time zero: \(S_0\).
- The option strike price: \(K\).
- The risk-free rate: \(rf\).
- The option maturity: \(T\).
- The volatility of the underlying stock \(\sigma\).
This is, \(BS = f(S_0, K, rf, T, \sigma)\). Here we set up the BS function and then evaluate it at \(BS = f(S_0=50, K=50, rf=0.05, T=1, \sigma=0.3)\).
# Black-Scholes call.
c.bs <- function(S, K, rf, TT, sigma) {
d1 <- (log(S / K) + (rf + sigma^2 / 2) * TT) / (sigma * sqrt(TT))
d2 <- d1 - sigma * sqrt(TT)
c <- S * pnorm(d1) - K * exp(-rf * TT) * pnorm(d2)
}
# Black-Scholes put.
p.bs <- function(S, K, rf, TT, sigma) {
d1 <- (log(S / K) + (rf + sigma^2 / 2) * TT) / (sigma * sqrt(TT))
d2 <- d1 - sigma * sqrt(TT)
p <- K * exp(-rf * TT) * pnorm(-d2) - S * pnorm(-d1)
}
# Calculate and report results.
c <- c.bs(50, 50, 0.05, 1, 0.3)
p <- p.bs(50, 50, 0.05, 1, 0.3)
data.frame("Black-Scholes call" = c, "Black-Scholes put" = p)## Black.Scholes.call Black.Scholes.put
## 1 7.115627 4.677099Now, let’s see how Black-Scholes option prices changes (call and put) when we leave \(S_0\) as a variable and the rest of the parameters as fixed: \(BS = f(S_0, K=50, rf=0.05, T=1, \sigma=0.3)\).
# Parameters.
K <- 50
rf <- 0.05
TT <- 1
sigma <- 0.3
S0 <- 50
# Sequence of S0.
S0.seq <- seq(from = 0, to = 100, length.out = 50)
# Evaluation of the c.bs and p.bs function.
c.S0 <- mapply(c.bs, S0.seq, K, rf, TT, sigma)
p.S0 <- mapply(p.bs, S0.seq, K, rf, TT, sigma)We produce as many c.S0 as values for S0.seq. Then, we can plot.
# See the results.
plot(S0.seq, c.S0, type = "l", lwd = 4, xlab = expression(paste(S[0])),
ylab = "Call option price, c")
Figure 2.1: Call price as a function of \(S_0\).
The call option price is an increasing function of \(S_0\). This makes sense because a right to buy the stock at a fix price in the future becomes more valuable as the current value of the stock increases. In other words, it is more expensive to lock the price of an expensive asset for a potential buyer.
See the case of the put option.
plot(S0.seq, p.S0, type = "l", lwd = 4, xlab = expression(paste(S[0])),
ylab = "Put option price, p")
Figure 2.2: Put price as a function of \(S_0\).
Here, it is cheaper to lock the price of an expensive asset for a potential seller. This makes sense because it looks unnecessary to lock the selling price of an expensive stock.
Let’s put both plots together.
# What is the interpretation of S0.star, and how did I found this value?
S0.star <- K * exp(-rf * TT)
c.S0.star <- mapply(c.bs, S0.star, K, rf, TT, sigma)
# Plot.
plot(S0.seq, c.S0, type = "l", lwd = 4, col = "blue", ylim = c(4, 8),
xlab = expression(paste(S[0])), ylab = "Option theoretical price",
xlim = c(40, 55))
lines(S0.seq, p.S0, lwd = 4, col = "red")
abline(v = S0.star, lty = 2)
abline(h = c.S0.star, lty = 2)
abline(h = c, lty = 2, col = "blue")
abline(h = p, lty = 2, col = "red")
abline(v = S0, lty = 2)
points(S0, c, pch = 19, col = "blue", cex = 3)
points(S0, p, pch = 19, col = "red", cex = 3)
points(S0.star, c.S0.star, pch = 19, cex = 3)
legend("topleft", legend = c("Call", "Put"),
col = c("blue", "red"), lwd = 3, bg = "white")
Figure 2.3: Call and put prices as a function of \(S_0\) (zoom).
Note that there is a single case in which both options (call and put) are the same. It is interesting to find out the value of \(S_0\) that makes \(c=p\). Clearly, this value is lower than \(S_0=50\), the point here is to find the mathematical expression of \(S_0\) that makes \(c=p\).
2.2 Option bounds.
Option values (or prices) are subject or constrained to upper and lower bounds. This basically mean that in theory option prices cannot exceed upper bounds and cannot decrease below lower bounds. If they do, then arbitrage opportunities arise.
We can extend our analysis by incorporating the corresponding call and put option bounds.
# Put lower bound.
lb.put <- pmax(K * exp(-rf * TT) - S0.seq, 0)
# Call lower bound.
lb.call <- pmax(S0.seq - K * exp(-rf * TT), 0)
# Again the S0.star. Where did this value comes from?
line.star <- S0.star * 2
# Plot.
par(pty = "s")
plot(S0.seq, c.S0, type = "l", lwd = 3, col = "blue", ylim = c(0, 100),
xlab = expression(paste(S[0])), ylab = "Call (blue), put (red)")
lines(S0.seq, p.S0, lwd = 3, col = "red")
abline(v = S0.star, lty = 2)
abline(v = 0, lty = 2)
abline(v = line.star, lty = 2, col = "green")
abline(h = line.star, lty = 2, col = "green")
# Put upper bound.
abline(h = K * exp(-rf * TT), lwd = 2, col = "red", lty = 2)
lines(S0.seq, lb.put, lwd = 2, col = "red", lty = 2)
# Call upper bound.
lines(S0.seq, S0.seq, lwd = 2, col = "blue", lty = 2)
lines(S0.seq, lb.call, lwd = 2, col = "blue", lty = 2)
Figure 2.4: Option values as a function of \(S_0\) including option bounds.
Note that the stock price range of values is from zero to positive infinity. However, this is not the case for the options. The option values depend on these bounds. We can illustrate the same plot in a cleaner way, to make emphasis on a geometric approach.
par(pty = "s")
plot(S0.seq, lb.put, type = "l", lwd = 2, col = "red", ylim = c(0, 100),
xlab = expression(paste(S[0])), ylab = "Call (blue), put (red)")
abline(v = S0.star, lty = 2)
abline(v = 0, lty = 2)
abline(v = line.star, lty = 2)
abline(h = line.star, lty = 2)
abline(h = S0.star, lwd = 2, col = "red")
lines(S0.seq, S0.seq, lwd = 2, col = "blue")
lines(S0.seq, lb.call, lwd = 2, col = "blue")
# Enumerate the areas.
text(S0.star / 2, (S0.star + (S0.star / 2)), "1")
text(60, 80, "2")
text(80, 60, "3")
text(10, 20, "4")
text(25, 35, "5")
text(25, 10, "6")
text(40, 20, "7")
text(60, 35, "8")
text(80, 20, "9")
Figure 2.5: Call and put option bounds: geometric inspection (or ‘the puzzle’).
Every number from 1 to 9 is assigned to a specific geometric shape in the plot above. This is, a square as in 1, and triangles of two different shapes in the rest of the cases (2 to 9). It is interesting to analyze whether each area represents a possible price range for call and/or put options.
Now, let’s see how Black-Scholes option prices changes (call and put) when we leave \(S_0\) as a variable, three levels of \(\sigma\) and the rest of the parameters as fixed: \(BS = f(S_0, K=50, rf=0.05, T=1, \sigma)\). We can show this in a 2 dimension plot, including the correspondent bounds.
# Evaluate the function.
c.S0.s1 <- mapply(c.bs, S0.seq, K, rf, TT, sigma)
c.S0.s5 <- mapply(c.bs, S0.seq, K, rf, TT, sigma * 5)
c.S0.s15 <- mapply(c.bs, S0.seq, K, rf, TT, sigma * 15)
# Plot.
plot(S0.seq, c.S0.s1, type = "l", lwd = 3, xlab = expression(paste(S[0])),
ylab = "Call option price, c")
lines(S0.seq, c.S0.s5, lwd = 3, col = "blue") # call with higher sigma
lines(S0.seq, c.S0.s15, lwd = 3, col = "red") # call with even higher sigma
lines(S0.seq, lb.call, lwd = 3, lty = 2) # lower bound
lines(S0.seq, S0.seq, lwd = 3, lty = 2) # upper bound
legend("topleft",legend=expression(sigma == 0.3, sigma == 1.5,
sigma == 4.5, "call bounds"),
text.width = 12, col = c("black", "blue", "red", "black"), lwd = 2,
lty = c(1, 1, 1, 2), bg = "white", cex = 0.8)
Figure 2.6: Call option changes as \(S_0\) and \(\sigma\) increases: the role of bounds.
And the same for the put option.
# Evaluate the function.
p.S0.s1 <- mapply(p.bs, S0.seq, K, rf, TT, sigma)
p.S0.s5 <- mapply(p.bs, S0.seq, K, rf, TT, sigma * 5)
p.S0.s15 <- mapply(p.bs, S0.seq, K, rf, TT, sigma * 15)
# Plot.
plot(S0.seq, p.S0.s1, type = "l", lwd = 3, xlab = expression(paste(S[0])),
ylab = "Pall option price, p")
lines(S0.seq, p.S0.s5, lwd = 3, col = "blue") # put with higher sigma
lines(S0.seq, p.S0.s15, lwd = 3, col = "red") # put with even higher sigma
lines(S0.seq, lb.put, lwd = 3, lty = 2) # lower bound
lines(S0.seq, rep(S0.star, 50) , lwd = 3, lty = 2) # upper bound
legend("bottomleft", legend=expression(sigma == 0.3, sigma == 1.5,
sigma == 4.5, "put bounds"),
text.width = 12, col = c("black", "blue", "red", "black"), lwd = 2,
lty = c(1, 1, 1, 2), bg = "white", cex = 0.8)
Figure 2.7: Put option changes as \(S_0\) and \(\sigma\) increase: the role of bounds.
# Evaluate the function.
p.S0.s1 <- mapply(p.bs, S0.seq, K, rf, TT, sigma)
p.S0.s5 <- mapply(p.bs, S0.seq, K, rf, TT, sigma * 5)
p.S0.s15 <- mapply(p.bs, S0.seq, K, rf, TT, sigma * 15)
# Plot.
plot(S0.seq, p.S0.s5, type = "l", lwd = 3, xlab = expression(paste(S[0])),
ylab = "Pall option price, p", col = "blue",
xlim = c(0, 1), ylim = c(47, 47.8))
lines(S0.seq, p.S0.s15, lwd = 3, col = "red") # put with even higher sigma
lines(S0.seq, lb.put, lwd = 1, lty = 2) # lower bound
lines(S0.seq, rep(S0.star, 50) , lwd = 1, lty = 2) # upper bound
Figure 2.8: Put option changes as \(S_0\) and \(\sigma\) increase: the role of bounds (zoom).
Then, regardless of the extreme value of parameters, option bounds represent the maximum and minimum values for the call and put options. These option prices are theoretical as they are calculated by implementing the Black-Scholes formula. We normally evaluate whether the option market prices violates these bounds because this will flag an opportunity to implement an arbitrage strategy to generate a risk-free profit.
A different way to illustrate the option properties is by looking at a plane. Here, we let \(S_0\) and \(K\) as a free variables and remain the rest as fixed: \(BS = f(S_0, K, rf=0.05, T=1, \sigma = 0.3)\).
# K as a variable.
K.seq <- S0.seq
# Create the empty matrix.
c.S0.s <- matrix(0, nrow = 50, ncol = 50)
# Fill the empty matrix.
for(i in 1:50){ # Is there an easier way to do this?
for(j in 1:50){
c.S0.s[i, j] <- c.bs(S0.seq[i], K.seq[j], rf, TT, sigma) } }
# Plot.
c.S0.s.plot <- persp(S0.seq, K.seq, c.S0.s, zlab = "Call",
xlab = "S_0",
ylab = "K", theta = -60, phi = 10, expand = 0.5, col = "orange",
shade = 0.2, ticktype = "detailed")
points(trans3d(S0, K, c, c.S0.s.plot), cex = 2, pch = 19, col = "blue")
Figure 2.9: Call value as \(S_0\) and \(K\) change, a plane view.
An alternative view is by showing a contour plot. A contour plot is a graphical technique for representing a 3-dimensional surface by plotting constant \(z\) slices (option prices), called contours, on a 2-dimensional format. That is, given a value for \(z\), lines are drawn for connecting the \((x, y)\) coordinates (\(S_0\), \(K\)) where that \(z\) value occurs.
contour(S0.seq, K.seq, c.S0.s, xlab = expression(paste(S[0])),
ylab = "K", lwd = 2, nlevels = 20)
points(S0, K, pch = 19, col = "blue", cex = 2)
Figure 2.10: Call value as \(S_0\) and \(K\) change, a contour view.
Note that the blue circle \((S_0=50, K=50)\) is between the contour line 5 and 10, this makes sense because the call option at \(call(S_0=50, K=50, rf=0.05, T=1, \sigma = 0.3) = 7.115627\).
plot_ly(type = "surface" , x = S0.seq, y = K.seq , z = c.S0.s ) %>%
layout(scene = list(xaxis = list(title = "Stock price"),
yaxis = list(title = "K",
zaxis = list(title = "Call price")))) %>%
hide_colorbar()Figure 2.11: Call value as \(S_0\) and \(K\) change, an interactive view.
An interesting case is the value of the put as a function of the time to maturity. Let’s see the simplest case.
TT.seq <- seq(from = 1, to = 25, length.out = 50)
p.TT <- mapply(p.bs, S0, K, rf, TT.seq, sigma)
plot(TT.seq, p.TT, type = "l", lwd = 4, ylab = "Put value", xlab = "T")
points(TT, p, col = "red", cex = 2, pch = 19)
Figure 2.12: European put as \(T\) increases: the mysterious hump.
The put option price first increases as the time to maturity increases and after reaching a maximum, the put option price decreases. What is the mathematical expression for the time to maturity which makes the put option price maximum? What is the reason or logic behind this mysterious hump? Those are interesting questions to address.
Let’s illustrate the case when time to maturity and the stock price at time zero changes.
p.S0.T <- matrix(0, nrow = 50, ncol = 50)
for(i in 1:50) { # Is there an easier way to do this?
for(j in 1:50) {
p.S0.T[i, j] <- p.bs(S0.seq[i], K, rf, TT.seq[j], sigma) } }
p.S0.T.plot <- persp(S0.seq, TT.seq, p.S0.T, zlab = "Put", xlab = "S0",
ylab = "T", theta = 330, phi = 10, expand = 0.5,
col = "orange", shade = 0.2, ticktype = "detailed")
points(trans3d(S0, TT, p, p.S0.T.plot), cex = 2, pch = 19, col = "blue")
Figure 2.13: European put value as \(S_0\) and \(T\) change: a plane.
That is not quite clear. So, here is the contour view.
contour(S0.seq, TT.seq, p.S0.T, xlab = expression(paste(S[0])),
ylab = "T", lwd = 2, nlevels = 40, drawlabels = TRUE)
points(S0, TT, pch = 19, col = "blue", cex = 2)
Figure 2.14: Put value as \(S_0\) and \(T\) change, a contour view, or the hypnotic plot.
Now the mysterious hump is clearer than before. Remember the value of \(put(S_0=50, K=50, rf=0.05, T=1, \sigma = 0.3) = 4.677099\).
plot_ly(type = "surface" , x = S0.seq, y = TT.seq , z = p.S0.T ) %>%
layout(scene = list(xaxis = list(title = "Stock price"),
yaxis = list(title = "T",
zaxis = list(title = "Put price")))) %>%
hide_colorbar()Figure 2.15: Put value as \(S_0\) and \(T\) change, an interactive view.
2.3 Put-call parity.
The main idea behind the put-call parity is to understand how call and put option prices are related as today, at \(t=0\). In Hull, the procedure to derive the put-call parity starts with the definition of two portfolios: (1) a call and a bond; (2) a put and a stock. Then, derive the corresponding payoff of each portfolio at maturity. Doing this is relatively easy because we do not need a valuation method as we know the payoff functions for options, bonds and stocks. Given that we can demonstrate that the value of two different portfolios are the same at \(T\), then we can conclude that these two portfolios are worth the same at \(t=0\) as well.
The put-call parity is: \(c+Ke^{-rT}=p+S_0\).
We can verify that this equation holds.
\(p=c+Ke^{-rT}-S_0 \rightarrow p=7.115627+50e^{-0.05\times 1}- 50 \rightarrow p=4.677098\).
\(c=p+S_0-Ke^{-rT} \rightarrow c=4.677098+50-50e^{-0.05\times 1} \rightarrow c=7.115627\).
Let’s take one step back and demonstrate that both portfolios are worth the same at maturity. We first define the assets payoffs.
ST.seq <- seq(from = 0, to = 150, length.out = 50)
cT <- pmax(ST.seq - K, 0) # call payoff.
pT <- pmax(K - ST.seq, 0) # put payoff.
BT <- rep(K, 50) # bond payoff (at maturity).
pc <- pmax(ST.seq, K) # this will be important.Then we plot.
par(pty = "s")
par(mfrow = c(1, 2), oma = c(0, 0, 2, 0))
par(pty = "s")
# Portfolio A.
plot(ST.seq, cT, type = "l", ylab = "Payoff",
xlab = expression(paste(S[T])), lwd = 2,
lty = 2, col = "blue", ylim = c(0, 150))
lines(ST.seq, (cT + BT), lwd = 6)
lines(ST.seq, BT, lwd = 2, lty = 2, col = "red")
legend("topleft", legend = c("Call option", "Zero coupon bond",
"Total (Portfolio A)"),
col = c("blue", "red", "black"), lwd = c(2, 2, 2), lty = c(2, 2, 1),
bg = "white", cex = 0.7)
par(pty = "s")
# Portfolio C.
plot(ST.seq, (pT + ST.seq), lwd = 6, type = "l", ylab = "",
xlab = expression(paste(S[T])), ylim = c(0, 150))
lines(ST.seq, ST.seq, lwd = 2, lty = 2, col = "orange")
lines(ST.seq, pT, lwd = 2, lty = 2, col = "purple")
legend("topleft", legend = c("Put option", "Share", "Total (Portfolio C)"),
col = c("purple", "orange", "black"), lwd = 2,
lty = c(2, 2, 1), bg = "white", cex = 0.7)
Figure 2.16: Portfolio A and C payoffs, the put-call parity.
The point here is that the black line (total) is the same in both cases. This is why we argue that the payoffs of both portfolios are worth the same. If this is so, then they are worth the same in time \(t=0\).
In sum, this is the payoff of portfolios A and C: \(\mathrm{max}(S_T,K)\). This is already stored in pc.
par(pty = "s")
plot(ST.seq, pc, type = "l", ylab = "Payoff",
xlab = expression(paste(S[T])), lwd = 5, ylim = c(0, 150))
Figure 2.17: Portfolio A and C payoffs, \(\mathrm{max}(S_T,K)\).
Nice. We can manipulate the put-call parity to create a synthetic stock \(S_0\). This is:
\(c+Ke^{-rT}=p+S_0 \rightarrow S_0 = c+Ke^{-rT}-p\).
Graphically:
par(pty = "s")
plot(ST.seq, (cT + BT - pT), type = "l", ylab = "Payoff",
xlab = expression(paste(S[T])),
lwd = 5, ylim = c(-50, 100), xlim = c(0, 150), col = "orange")
lines(ST.seq, cT, lwd = 2, lty = 2, col = "blue") # positive call.
lines(ST.seq, BT, lwd = 2, lty = 2, col = "red") # positive bond.
lines(ST.seq, -pT, lwd = 2, lty = 2, col = "purple") # negative put.
legend("bottomright", legend = c("Call option", "Zero coupon bond",
"Put option", "Total (stock)"),
col = c("blue", "red", "purple", "orange"), lwd = 2,
lty = c(2, 2, 2, 1), bg = "white", cex = 0.65)
Figure 2.18: A synthetic stock \(S_0\).
So, we created a stock that did not exist with a call, a bond and a put. We can manipulate the put-call parity to create a synthetic bond \(Ke^{-rT}\). This is:
\(c+Ke^{-rT}=p+S_0 \rightarrow Ke^{-rT} = p+S_0-c\).
Graphically:
par(pty = "s")
plot(ST.seq, (-cT + ST.seq + pT), type = "l", ylab = "Payoff",
xlab = expression(paste(S[T])),
lwd = 5, ylim = c(-50, 100), xlim = c(0, 150), col = "orange")
lines(ST.seq, -cT, lwd = 2, lty = 2, col = "blue")
lines(ST.seq, ST.seq, lwd = 2, lty = 2, col = "red")
lines(ST.seq, pT, lwd = 2, lty = 2, col = "purple")
legend("bottomleft", legend = c("Call option (short)", "Stock",
"Put option", "Total (bond)"),
col = c("blue", "red", "purple", "orange"), lwd = 2,
lty = c(2, 2, 2, 1), bg = "white", cex = 0.65)
Figure 2.19: A synthetic bond \(Ke^{-rT}\).
There is an alternative to the Black-Scholes formula we indirectly review before. Next section introduces the binomial trees.
3 Binomial trees.
Binomial trees are a flexible valuation method for options because we can use them to value not only European but also American options. This section is based in Hull (2015) chapter 13.
3.1 Implementation.
The function is the following.
bin <- function(S0, K, sigma, TM, r, steps) {
# the parameters
dt <- TM / steps
u <- exp(sigma * sqrt(dt))
d <- exp(-sigma * sqrt(dt))
a <- exp(r * dt)
p <- (a - d) / (u - d)
S <- matrix(0, steps + 1, steps)
pam <- S
peu <- S
cam <- S
ceu <- S
# the stock price process
for(i in 1:steps) {
j <- i + 1
do <- seq(0, i)
up <- rev(do)
S[(1:j), i] <- S0 * (u ^ up) * (d ^ do) }
# the option prices at maturity.
peu[(1:(steps + 1)), steps] <- pmax(K - S[, i], 0)
pam[(1:(steps + 1)), steps] <- pmax(K - S[, i], 0)
ceu[(1:(steps + 1)), steps] <- pmax(S[, i] - K, 0)
cam[(1:(steps + 1)), steps] <- pmax(S[, i] - K, 0)
# the binomial method to price stock options.
for(j in steps:1) { # this is a reverse loop from steps to 1.
cd <- (seq(steps:1)) # every round we compute less option prices.
for(i in 1:cd[j]) { # option prices per step.
peu[i, (j - 1)] <- exp(-r * dt) * (p * peu[i, j] +
(1 - p) * peu[(i + 1), j])
ceu[i, (j - 1)] <- exp(-r * dt) * (p * ceu[i, j] +
(1 - p) * ceu[(i + 1), j])
pam[i, (j - 1)] <- max((K - S[i, (j - 1)]), exp(-r * dt) *
(p * pam[i, j] + (1 - p) * pam[(i + 1), j]))
cam[i, (j - 1)] <- max((S[i, (j - 1)] - K), exp(-r * dt) *
(p * cam[i, j] + (1 - p) * cam[(i + 1), j])) } }
# This is the final step in the binomial tree.
p.eu <- exp(-r * dt) * (p * peu[1, 1] + (1 - p) * peu[2, 1])
c.eu <- exp(-r * dt) * (p * ceu[1, 1] + (1 - p) * ceu[2, 1])
p.am <- exp(-r * dt) * (p * pam[1, 1] + (1 - p) * pam[2, 1])
c.am <- exp(-r * dt) * (p * cam[1, 1] + (1 - p) * cam[2, 1])
# Results.
option <- data.frame(c.eu, p.eu, c.am, p.am)
option
}We can evaluate the function above to see how the price change depending on the number of steps in the binomial tree. Let’s evaluate the bin() function assuming \(S_0=50\), \(K=50\), \(rf=0.05\), \(T=1\), and \(\sigma=0.3\) at different time steps.
b1t <- bin(S0, K, sigma, TT, rf, 1)
b4t <- bin(S0, K, sigma, TT, rf, 4)
b20t <- bin(S0, K, sigma, TT, rf, 20)
b50t <- bin(S0, K, sigma, TT, rf, 50)
b200t <- bin(S0, K, sigma, TT, rf, 200)
b500t <- bin(S0, K, sigma, TT, rf, 500)See the results.
Black.Scholes <- data.frame(c.eu = c, p.eu = p, c.am = NA, p.am = NA)
bin.bs <- rbind("Binomial (1 step)" = b1t, "Binomial (4 steps)" = b4t,
"Binomial (20 steps)" = b20t, "Binomial (50 steps)" = b50t,
"Binomial (200 steps)" = b200t,
"Binomial (500 steps)" = b500t,
"Black-Scholes" = Black.Scholes)
bin.bs## c.eu p.eu c.am p.am
## Binomial (1 step) 8.481986 6.043457 8.481986 6.043457
## Binomial (4 steps) 6.762001 4.323472 6.762001 4.767526
## Binomial (20 steps) 7.042462 4.603934 7.042462 4.898985
## Binomial (50 steps) 7.086241 4.647713 7.086241 4.921038
## Binomial (200 steps) 7.108267 4.669738 7.108267 4.931581
## Binomial (500 steps) 7.112682 4.674153 7.112682 4.933664
## Black-Scholes 7.115627 4.677099 NA NAIn the extreme, as the number of steps increases, the binomial method converges to the Black-Scholes method. Let’s explore these differences.
b1 <- mapply(bin, S0.seq, K, sigma, TT, rf, 1)
b2 <- mapply(bin, S0.seq, K, sigma, TT, rf, 2)
b3 <- mapply(bin, S0.seq, K, sigma, TT, rf, 3)
b4 <- mapply(bin, S0.seq, K, sigma, TT, rf, 4)
b20 <- mapply(bin, S0.seq, K, sigma, TT, rf, 20)Here, we compare the value of the binomial method and the Black-Scholes at different time steps as a function of \(S_0\).
plot(S0.seq, b1[1, ], type = "l", col = "red", lwd = 5,
xlab = expression(paste(S[0])), ylab = "European call")
lines(S0.seq, c.S0.s1, col = "black", lwd = 5)
legend("topleft", legend = c("1 step binomial", "Black-Scholes"),
col = c("red", "black"), lwd = 2, bg = "white")
Figure 3.1: One-step binomial and Black-Scholes comparison.
There are some differences between the binomial and the Black-Scholes. Let’s zoom to see the differences clearer.
plot(S0.seq, b1[1,], type = "l", ylim = c(0, 4), xlim = c(30, 45),
col = "green", lwd = 2, xlab = expression(paste(S[0])),
ylab = "European call")
lines(S0.seq, b2[1,], col = "purple", lwd = 2)
lines(S0.seq, b3[1,], col = "orange", lwd = 2)
lines(S0.seq, b4[1,], col = "red", lwd = 2)
lines(S0.seq, b20[1,], col = "black", lwd = 2)
lines(S0.seq, c.S0.s1, col = "black", lwd = 2)
legend("topleft", legend = c("1 step binomial", "2 steps binomial",
"3 steps binomial", "4 steps binomial",
"20 steps binomial", "Black-Scholes"),
col = c("green", "purple", "orange", "red", "black", "black"),
lwd = 2, bg = "white")
Figure 3.2: Binomial and Black-Scholes convergence, a zoom view.
Not clear yet. Let’s see a panel view.
par(mfrow = c(2, 2), mai = c(0.7, 0.4, 0.4, 0.4))
par(pty = "s")
plot(S0.seq, c.S0.s1, type = "l", ylim = c(0, 4), xlim = c(30, 45),
col = "black", lwd = 6, xlab = "", ylab = "European call")
lines(S0.seq, b2[1,], col = "purple", lwd = 2)
legend("topleft", legend = c("2-steps"),
bg = "white", bty = "n")
par(pty = "s")
plot(S0.seq, c.S0.s1, type = "l", ylim = c(0, 4), xlim = c(30, 45),
col = "black", lwd = 6, xlab = "", ylab = "")
lines(S0.seq, b3[1,], col = "orange", lwd = 2)
legend("topleft", legend = c("3-steps"),
bg = "white", bty = "n")
par(pty = "s")
plot(S0.seq, c.S0.s1, type = "l", ylim = c(0, 4), xlim = c(30, 45),
col = "black", lwd = 6, xlab = expression(paste(S[0])),
ylab = "European call")
lines(S0.seq, b4[1,], col = "red", lwd = 2)
legend("topleft", legend = c("4-steps"), bg = "white", bty = "n")
par(pty = "s")
plot(S0.seq, c.S0.s1, type = "l", ylim = c(0, 4), xlim = c(30, 45),
col = "black", lwd = 6, xlab = expression(paste(S[0])), ylab = "")
lines(S0.seq, b20[1,], col = "grey", lwd = 2)
legend("topleft", legend = c("20-steps"),
bg = "white", bty = "n")
Figure 3.3: Binomial and Black-Scholes convergence, a panel view.
As stated above, the binomial method converges to the Black-Scholes as time steps increase. We can also create a function to visualize the price path of the stock from \(S_0\) to \(S_T\) given the assumptions of the binomial method.
3.2 Stock price paths.
In order to value option prices, we first need to understand the evolution of the underlying (in this case the stock price) from time zero to maturity, this is from \(S_0\) to \(S_T\). The binomial method assumes that the price can increase or decrease with a certain probability in each time step. Here, we can show how the binomial tree method assumes this stock price evolution.
Let’s assume 10-steps, from \(t=0.1, t=0.2,..., t=1\)
# Function to generate stock prices paths given the binomial method.
S.paths <- function(S0, sigma, TM, steps) {
dt <- TM / steps
u <- exp(sigma * dt^0.5) # Here we set u and d as a function of sigma.
d <- exp(-sigma * dt^0.5)
S <- matrix(0, (steps + 1), (steps + 1))
S[1, 1] <- S0
for(i in 2:(steps + 1)) {
do = seq(0, i - 1)
up = rev(do) # rev provides a reversed version of its argument.
S[(1:i), i] = S0 * (u ^ up) * (d ^ do) }
S }
# Evaluate the function.
Spaths <- S.paths(50, 0.3, 1, 10)
# A table.
colnames(Spaths) <- c(0, cumsum(rep(1/10, 10)))
round(Spaths, 2)## 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
## [1,] 50 54.98 60.45 66.46 73.08 80.35 88.34 97.13 106.80 117.43 129.12
## [2,] 0 45.47 50.00 54.98 60.45 66.46 73.08 80.35 88.34 97.13 106.80
## [3,] 0 0.00 41.36 45.47 50.00 54.98 60.45 66.46 73.08 80.35 88.34
## [4,] 0 0.00 0.00 37.62 41.36 45.47 50.00 54.98 60.45 66.46 73.08
## [5,] 0 0.00 0.00 0.00 34.21 37.62 41.36 45.47 50.00 54.98 60.45
## [6,] 0 0.00 0.00 0.00 0.00 31.11 34.21 37.62 41.36 45.47 50.00
## [7,] 0 0.00 0.00 0.00 0.00 0.00 28.30 31.11 34.21 37.62 41.36
## [8,] 0 0.00 0.00 0.00 0.00 0.00 0.00 25.74 28.30 31.11 34.21
## [9,] 0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 23.41 25.74 28.30
## [10,] 0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 21.29 23.41
## [11,] 0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 19.36This does not look like a typical binomial tree. In fact, it is not very clear whether a given price corresponds to an increase or decrease from a previous time step. We can make a few arrangements to visualize this as a tree.
S.paths <- function(S0, sigma, TM, steps) {
dt <- TM / steps
u <- exp(sigma * dt^0.5)
d <- exp(-sigma * dt^0.5)
S <- matrix(0, 2 * ((steps + 1)), (steps + 1))
S2 <- matrix(NA, ((2 * steps) + 2), (steps + 1))
S2[(steps + 1), 1] <- S0
for(i in 2:(steps + 1)) {
do = seq(0, i - 1)
up = rev(do) # rev provides a reversed version of its argument.
S[(1:i), i] = S0 * (u ^ up) * (d ^ do)
x = rep(NA, i) # These are the NA between stock prices.
r = rev(c(seq(0, (steps - 1)), 0)) # These creates the blank spaces.
# Here we combine NA and stock prices for each column.
S2[(1 + r[i]):((2 * i) + r[i]), i] = as.numeric(rbind(S[(1:i), i], x))
}
S2 }
# Evaluate the function.
Spaths <- S.paths(50, 0.3, 1, 10)
# A table.
colnames(Spaths) <- round(c(0, cumsum(rep(1/10, 10))), 2)
round(Spaths, 2)## 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
## [1,] NA NA NA NA NA NA NA NA NA NA 129.12
## [2,] NA NA NA NA NA NA NA NA NA 117.43 NA
## [3,] NA NA NA NA NA NA NA NA 106.80 NA 106.80
## [4,] NA NA NA NA NA NA NA 97.13 NA 97.13 NA
## [5,] NA NA NA NA NA NA 88.34 NA 88.34 NA 88.34
## [6,] NA NA NA NA NA 80.35 NA 80.35 NA 80.35 NA
## [7,] NA NA NA NA 73.08 NA 73.08 NA 73.08 NA 73.08
## [8,] NA NA NA 66.46 NA 66.46 NA 66.46 NA 66.46 NA
## [9,] NA NA 60.45 NA 60.45 NA 60.45 NA 60.45 NA 60.45
## [10,] NA 54.98 NA 54.98 NA 54.98 NA 54.98 NA 54.98 NA
## [11,] 50 NA 50.00 NA 50.00 NA 50.00 NA 50.00 NA 50.00
## [12,] NA 45.47 NA 45.47 NA 45.47 NA 45.47 NA 45.47 NA
## [13,] NA NA 41.36 NA 41.36 NA 41.36 NA 41.36 NA 41.36
## [14,] NA NA NA 37.62 NA 37.62 NA 37.62 NA 37.62 NA
## [15,] NA NA NA NA 34.21 NA 34.21 NA 34.21 NA 34.21
## [16,] NA NA NA NA NA 31.11 NA 31.11 NA 31.11 NA
## [17,] NA NA NA NA NA NA 28.30 NA 28.30 NA 28.30
## [18,] NA NA NA NA NA NA NA 25.74 NA 25.74 NA
## [19,] NA NA NA NA NA NA NA NA 23.41 NA 23.41
## [20,] NA NA NA NA NA NA NA NA NA 21.29 NA
## [21,] NA NA NA NA NA NA NA NA NA NA 19.36
## [22,] NA NA NA NA NA NA NA NA NA NA NAGraphically:
nona <- ! is.na(Spaths)
bin_paths <- data.frame(st = rep(1:11, 1:11)/10-0.1, pa = Spaths[nona])
plot(bin_paths$st, bin_paths$pa, pch = 19, ylab = expression(paste(S[t])),
xlab = "t")
Figure 3.4: \(S_t\) binomial paths from \(S_0\) to \(S_T\).
Nice. Evaluate the function for a single 2-step case: \(S_0=50\), \(K=52\), \(\sigma = 0.3\), \(T=2\), \(rf=0.05\).
# bin <- function(S0, K, sigma, TM, r, steps) {
b2 <- mapply(bin, 50, 52, 0.3, 2, 0.05, 2)
b2## [,1]
## c.eu 9.194163
## p.eu 6.245708
## c.am 9.194163
## p.am 7.428402There is no difference between American and European call options theoretical prices as there are no incentives to exercise American call options early. The case of put options is different as incentives to exercise American put options early may arise. As a consequence, American put options are in general more expensive than European put options.
Let’s evaluate this difference in a 50-step case: \(K=52\), \(\sigma = 0.3\), \(T=2\), \(rf=0.05\), and different values of \(S_0\).
S.seq <- seq(0, 80, 0.5)
AEoptions <- mapply(bin, S.seq, 52, 0.3, 1, 0.05, 50)
plot(S.seq, unlist(AEoptions[4,]) - unlist(AEoptions[2,]),
type = "l", lwd = 5, main = "The value of being American.",
xlab = expression(paste(S[0])),
ylab = "American minus European put")
abline(v = 52, lty = 2)
abline(h = 0, lty = 2)
Figure 3.5: The value of being American; or the cost of being European.
3.3 Real vs. risk-neutral world.
This section is based on Hull (2015) section 13.1. A stock price is currently $20, and it is known that at the end of 3 months it will be either $22 or $18. We are interested in valuing a European call option to buy the stock for $21 in 3 months. The risk-neutral probability \(p\) for the case of \(rf=0.12\), \(T=3/12\), \(u=1.1\) and \(d=0.9\):
\(p = \frac{e^{rT}-d}{u-d} \rightarrow p = \frac{e^{0.12\times(3/12)}-0.9}{1.1-0.9} \rightarrow p = 0.6522727\).
# This is the substitution of equation 13.3.
rnp <- (exp(0.12*(3/12))-0.9)/(1.1-0.9)
rnp## [1] 0.6522727Following Hull (2015), this option will have one of two values at the end of the 3 months. If the stock price turns out to be $22, the value of the option will be $1; if the stock price turns out to be $18, the value of the option will be zero. This is the European call option expected value at time \(T\):
\([p\times(\$1)] + [(1-p)\times \$0]\)
\(\rightarrow [0.6522727\times(\$1)] + [(1-0.6522727)\times \$0] = \$0.6522727\).
option_T <- rnp * 1 + (1 - rnp) * 0
option_T## [1] 0.6522727We can calculate the European call option value at \(t=0\) according to the risk-neutral valuation. This option value \(f\) is valid in the real world, not only in the risk-neutral world.
\(f = \$0.6522727 \times e^{-rT} \rightarrow f = \$0.6522727 \times e^{-0.12\times3/12} \rightarrow f = \$0.6329951\).
option_0 <- option_T * exp(-0.12*(3/12))
option_0## [1] 0.6329951Suppose that, in the real world, we know the expected return on the stock is 16%. This is easily calculated as we have historical information of stock returns. It is also a reasonable assumption as this is higher than the risk-free rate. Then, the probability of an up movement in the real world \(p^*\) is:
\(p^* = \frac{e^{rT}-d}{u-d} \rightarrow p^* = \frac{e^{0.16\times(3/12)}-0.9}{1.1-0.9} \rightarrow p^* = 0.7040539\).
# Real world probability.
p.real <- (exp(0.16*(3/12))-0.9)/(1.1-0.9)
p.real## [1] 0.7040539This is the option expected value in the real world at time \(T\):
\([p^*\times(\$1)] + [(1-p^*)\times \$0]\)
\(\rightarrow 0.7040539\times(\$1) + (1-0.7040539)\times \$0 = \$0.7040539\).
option_T_real <- p.real*1 + (1-p.real)*0
option_T_real## [1] 0.7040539Can we calculate the present value of $0.7040539? We only need the option discount factor in the real world. The problem is that we do not know this rate, it should be higher than the risk-free rate, even higher than the the expected return on the stock. It is not easy to know the correct discount rate of the option to apply to the expected payoff in the real world. Using risk-neutral valuation solves this problem because we know that in a risk-neutral world the expected return on all assets (and therefore the discount rate to use for all expected payoffs) is the risk-free rate.
What can we do to find out the option discount factor? Since we know the correct value of the option is \(f=\$0.6329951\) in both real and risk-neutral world, we can deduce the real-world discount rate \(r\) of this European call option.
We know this is true: \(\$0.6329951 = \$0.704053e^{-r\times3/12}\).
Then, solve for \(r\): \(\rightarrow \frac{\$0.6329951}{\$0.704053} = e^{-r\times3/12}\),
\(\rightarrow log{\frac{\$0.6329951}{\$0.704053}} = -r\times3/12\),
\(\rightarrow r = -log{\frac{\$0.6329951}{\$0.704053}}\times12/3\),
\(\rightarrow r=0.4255688\).
option.df <- -log(option_0/option_T_real)*(12/3)
option.df## [1] 0.4255688The implied real-world discount rate for the option is 42.55688%. Let’s put everything in a function.
option.df.fun <- function(r) {
p.real <- (exp(r*(3/12))-0.9)/(1.1-0.9)
option_T_real <- p.real*1 + (1-p.real)*0
option.df <- -log(option_0/option_T_real)*(12/3)
option.df
}Evaluate for known values:
option.df.fun(0.12)## [1] 0.12option.df.fun(0.16)## [1] 0.4255688Everything looks correct.
# Evaluate the function in a range of values.
r.seq <- seq(0.11, 0.17, 0.001)
df <- mapply(option.df.fun, r.seq)
plot(r.seq, df, type = "l", lwd = 5,
xlab = "Expected return of the stock (the underlying)",
ylab = "Expected return of the option")
points(0.12, 0.12, pch = 19, col = "blue", cex = 3)
points(0.16, 0.4255688, pch = 19, col = "red", cex = 3)
abline(v = 0.12, lty = 2, col = "blue", lwd = 2)
abline(h = option.df.fun(0.12), lty = 2, col = "blue", lwd = 2)
abline(v = 0.16, lty = 2, col = "red", lwd = 2)
abline(h = option.df.fun(0.16), lty = 2, col = "red", lwd = 2)
Figure 3.6: Red pill, blue pill.
As in the movie, The Matrix, the blue pill describes living life without knowing its meaning or running away from the truth in order to stay as is. This is equivalent to the risk-neutral world. The red pill is described as the solution for knowing the real truth in life. This is equivalent to the real-world.
3.4 Parrondo’s paradox.
The Parrondo’s paradox is controversial in the literature, we are not going to elaborate on that. However, it is indeed interesting as it looks counterintuitive: Can we combine two losing investments into a winner?
This video was created by Professor Humberto Barreto from DePauw University, Indiana, based on a now-defunct app from Alexander Bogomolny of the well-known maths site Cut the Knot. It illustrates the general idea of Parrondo’s paradox.
embed_url("https://youtu.be/cEuyfD2qVgQ")Consider asset A which has a known 6% gross return. In the context of a binomial tree: \([p\times(1+u)]+[(1-p)\times(1-d)]=1.06\). Let’s assume \(p=0.5\), so we have:
\(0.5(1+u)+0.5(1-d)=1.06\).
Asset A has a known 40% volatility, so:
\(0.4=\sqrt{0.5(1+u)^2+0.5(1-d)^2-1.06^2}\).
We have two equations and two unknown variables, so we can solve for \((1+u)\) and \((1-d)\): \((1+u)=1.46\) and \((1-d)=0.66\). In a period of 5 years, we would have a random path of ups (1.46) and downs (0.66) for asset A.
u = 0.46
d = 0.34
set.seed(2, sample.kind = "Rounding")
a <- sample(c(1+u, 1-d), 5, replace = TRUE)
a## [1] 1.46 0.66 0.66 1.46 0.66Graphically:
barplot(c(0, a-1), ylim = c(-0.6, 0.6), ylab = c("Returns: u or d."),
col = ifelse(c(0, a-1) < 0, "red", "blue"))
abline(h = 0)
Figure 3.7: 5-year simulation of \(u\) and \(d\) for asset A.
The evolution of $1 dollar invested in this 5 year period would look like this:
cumprod(a)## [1] 1.4600000 0.9636000 0.6359760 0.9285250 0.6128265Graphically:
plot(c(0:5), c(1, cumprod(a)), type = "b", lwd = 5,
ylab = c("Cumulative $ return"), xlab = c("Year"))
Figure 3.8: 5-year simulation of cumulative returns for asset A.
Even simpler and quicker, a \(\$1\) dollar invested at \(t=0\) would lead to \(\$0.6128265\) at \(t=5\).
prod(a)## [1] 0.6128265Let’s use a function now.
binomial_tree <- function(u, d, n) {
prod(sample(c(1+u, 1-d), n, replace = TRUE))
}See if it works. Calculate the asset A cumulative return of a 5-year and 30-year investment.
set.seed(2, sample.kind = "Rounding")
binomial_tree(u = 0.46, d = 0.34, n = 5)## [1] 0.6128265set.seed(2, sample.kind = "Rounding")
binomial_tree(u = 0.46, d = 0.34, n = 30)## [1] 2.805884It works, a \(\$1\) dollar invested at \(t=0\) would lead to \(\$0.6128265\) in 5 years as we show before. And now we know that a \(\$1\) dollar invested at \(t=0\) would lead \(\$2.805884\) in 30 years. Let’s view the 30 year case:
u = 0.46
d = 0.34
n = 30
set.seed(2, sample.kind = "Rounding")
a.30 <- c(1, cumprod(sample(c(1+u, 1-d), n, replace = TRUE)))
a.30 <- as.data.frame(cbind(year = c(0:30), ret = a.30))
ggplot(a.30, aes(x = year, y = ret)) +
geom_line(size = 1) +
geom_point(aes(x = 5, y = 0.6128265), size = 6) +
labs(y = "Cumulative return ($)", x = "Year",
subtitle = "Evolution of $1 invested in year 0.")
Figure 3.9: 30-year cumulative return asset A.
These results represent one single path. So, if we are interested in the most likely value of our investment we need to simulate many paths and then estimate the median. Let’s simulate 10,000 paths and see the most likely value of our \(\$1\) dollar investment in 30 years.
u = 0.46
d = 0.34
n = 30
set.seed(2, sample.kind = "Rounding")
a.30x10k <- replicate(10000,
cumprod(sample(c(1+u, 1-d), n, replace = TRUE))) |>
as.data.frame() |>
mutate(year = c(1:30)) |>
gather(V1:V10000, key = name, value = c.ret)See the results.
a.30x10k |>
ggplot(aes(x = year, y = c.ret, color = name)) +
geom_line(size = 1) +
geom_point(aes(x = 0, y = 1), col = "black", size = 4) +
theme(legend.position = "none", legend.title = element_blank()) +
labs(y = "Cumulative return ($)", x = "Year",
subtitle = "10,000 paths of the evolution of $1 invested in year 0.")
Figure 3.10: 30-year cumulative return asset A, many paths.
Most of these paths are concentrated in low cumulative returns by the end of the 30 years. In fact, 1,586 paths end up with a cumulative return lower than $1 by the end of the 30 years:
a.30x10k |>
group_by(name) |>
filter(c.ret < 1) |>
filter(length(year) == 30) |>
distinct(name) |>
nrow()## [1] 1586Let’s plot these 1,586 paths.
a.30x10k |>
group_by(name) |>
filter(c.ret < 1) |>
filter(length(year) == 30) |>
ggplot(aes(x = year, y = c.ret, color = name)) +
geom_line(size = 1) +
geom_point(aes(x = 0, y = 1), col = "black", size = 4) +
theme(legend.position = "none", legend.title = element_blank()) +
labs(y = "Cumulative return ($)", x = "Year",
subtitle = "1,586 paths of the evolution of $1 invested in year 0.")
Figure 3.11: 30-year cumulative return asset A, 1,586 selected paths.
The figure above looks still unclear as the concentration around $0 is still high. It is better to show the results in logarithm form.
a.30x10k |>
group_by(name) |>
filter(c.ret < 1) |>
filter(length(year) == 30) |>
ggplot(aes(x = year, y = log(c.ret), color = name)) +
geom_line(size = 1) +
geom_point(col = "black", alpha = 0.01, size = 4) +
geom_point(aes(x = 0, y = log(1)), col = "black", size = 4) +
theme(legend.position = "none",
legend.title = element_blank()) +
labs(y = "Log cumulative return", x = "Year",
subtitle = "Losing paths.")
Figure 3.12: 30-year log cumulative return asset A, 1,586 selected paths.
The figure now resembles a binomial tree. The black circles are more transparent for less concentrated values and less transparent for more concentrated values.
Let’s do a similar plot for all the 10,000 paths.
a.30x10k |>
ggplot(aes(x = year, y = log(c.ret), color = name)) +
geom_line(size = 1) +
geom_point(aes(x = 0, y = log(1)), size = 3, col = "black") +
geom_point(col = "black", alpha = 0.01, size = 3) +
theme(legend.position = "none", legend.title = element_blank()) +
labs(y = "Log cumulative returns", x = "Year",
subtitle = "10,000 paths of the evolution of $1 invested in year 0.")
Figure 3.13: 30-year Log cumulative return asset A, many paths.
Here it is easier to see that the paths actually follow a binomial structure. Now, let’s see the summary statistics for the 10,000 possible investment value at year 30.
set.seed(2, sample.kind = "Rounding")
summary(replicate(10000, binomial_tree(u = 0.46, d = 0.34, 30)))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.117 0.573 5.816 2.806 3559.180Not good news as our 30 year investment of \(\$1\) in this asset A with 6% gross return and 40% volatility leads to \(\$0.5733923\). Even if we drop the set.seed(), we end with the same result.
median(replicate(10000, binomial_tree(u = 0.46, d = 0.34, 30)))## [1] 0.5733923It would be a mistake to consider the mean $5.816 as the estimated value of our investment in asset A since the maximum value is very high and not likely to happen. Note that the mean is very close to \((\$1+0.06)^{30}=\$5.743491\).
(1+0.06)^30## [1] 5.743491The asset A is then a loosing investment, \(\$1\) would most likely lead to \(\$0.5733923\) in 30 years.
Now, let’s consider an asset B. Asset B has a known excess return of \(-0.1\%\). In the context of a binomial tree:
\(0.5(1+u)+0.5(1-d)=0.999\).
Asset B has no volatility, imagine is a Treasury Bill or a similar risk-free asset. Solving for \((1+u)\) and \((1-d)\) leads to \((1+u)=0.999\) and \((1-d)=0.999\). Asset B is for sure a loosing investment as well.
In a period of 5 years, we would have a random path of ups and downs in the context of a binomial path. Although in this case the result is the same as we have no risk.
b <- sample(c(1+(-0.001), 1-(0.001)), 5, replace = TRUE)
b## [1] 0.999 0.999 0.999 0.999 0.999The evolution of $1 dollar invested in this 5-year period would look like this:
cumprod(b)## [1] 0.999000 0.998001 0.997003 0.996006 0.995010Or simply:
prod(b)## [1] 0.99501As stated earlier, this asset B is a loosing investment. Let’s verify this for 5 and 30 years as we did before for the case of asset A:
binomial_tree <- function(u, d, n) {
prod(sample(c(1+u, 1-d), n, replace = TRUE))
}
set.seed(2, sample.kind = "Rounding")
binomial_tree(u = (-0.001), d = (0.001), 5)## [1] 0.99501binomial_tree(u = (-0.001), d = (0.001), 30)## [1] 0.970431It works, a \(\$1\) dollar invested at \(t=0\) would lead to \(\$0.99501\) in 5 years. And now we know that a \(\$1\) dollar invested at \(t=0\) would lead to \(\$0.970431\) in 30 years. These results represent one single path, but we do not need more as there is no risk. Let’s confirm:
set.seed(2, sample.kind = "Rounding")
median(replicate(10000, binomial_tree(u = -0.001, d = 0.001, 30)))## [1] 0.970431The result is the same. Asset B is a loosing investment as well. But, what if we combine both assets A and B in an equally weighted portfolio C?
set.seed(2, sample.kind = "Rounding")
median(replicate(1, binomial_tree(u = (0.46-0.001)/2,
d = (0.34+0.001)/2, 30)))## [1] 2.951197This is how we combine two losing investments into a winner.
The new asset C or portfolio C has a return of: \(0.06/2 - 0.001/2=0.0295%\) or \(2.95\%\). The volatility of C is \(\sqrt{0.5^2(0.4)^2+0.5^2(0)^2-0}=0.2\) or \(20\%\). Therefore, portfolio C has a lower return and risk than asset A. Given that A and B are uncorrelated, the diversification gain is high so we can combine two losing investments into a winner.
Let’s visualize how these assets behave. I have to set the seed to have nice results. But as we show before this works well when simulating many paths and computing the median to evaluate the expected cumulative return of each asset.
u = 0.46
d = 0.34
xx = 8
set.seed(xx, sample.kind = "Rounding")
a30 <- sample(c(1+u, 1-d), 30, replace = TRUE)
u = - 0.001
d = 0.001
b30 <- sample(c(1+u, 1-d), 30, replace = TRUE)
u = (0.46 - 0.001)/2
d = (0.34 + 0.001)/2
set.seed(xx, sample.kind = "Rounding")
c30 <- sample(c(1+u, 1-d), 30, replace = TRUE)Below, it is clear that asset C wins less when asset A wins. Similarly, asset C looses less when asset A losses. This is because asset C is less risky than asset A.
abc <- as.data.frame(cbind(
year = c(1:30), A = cumprod(a30),
B = cumprod(b30), C = cumprod(c30)))
head(abc)## year A B C
## 1 1 1.4600000 0.999000 1.229500
## 2 2 2.1316000 0.998001 1.511670
## 3 3 1.4068560 0.997003 1.253930
## 4 4 0.9285250 0.996006 1.040135
## 5 5 1.3556464 0.995010 1.278846
## 6 6 0.8947267 0.994015 1.060803abc |>
gather(A:C, key = name, value = c.ret) |>
ggplot(aes(x = year, y = c.ret, color = name)) +
geom_line(size = 1) +
geom_point(aes(x = 0, y = 1), size = 5, col = "black") +
labs(y = "Cumulative returns", x = "Year")
Figure 3.14: 30-year cumulative returns asset A, B and portfolio C.
set.seed(2, sample.kind = "Rounding")
bx <- (replicate(10000, binomial_tree(u = -0.001, d = 0.001, 30)))
set.seed(2, sample.kind = "Rounding")
cx <- (replicate(10000, binomial_tree(u = (0.46-0.001)/2,
d = (0.34+0.001)/2, 30)))
set.seed(2, sample.kind = "Rounding")
ax <-(replicate(10000, binomial_tree(u = 0.46, d = 0.34, 30)))
abcx <- as.data.frame(cbind(year = c(1:30), A = ax, B = bx, C = cx))
abcx |>
gather(A:C, key = name, value = c.ret) |>
ggplot(aes(x = name, y = log(c.ret), color = name)) +
geom_boxplot() +
theme(legend.position = "none") +
labs(y = "Log cumulative returns", x = "")
Figure 3.15: 30-year Log10 cumulative returns asset A, B and portfolio C, many paths.
Stock price paths are useful to value option stocks. Let’s review some stochastic processes in the next section.
library(knitr)
library(kableExtra)
library(plotly)4 Wiener processes.
Here we review stochastic process that are useful in finance, especially in the Black-Scholes context.
4.1 The basic process.
A simple Wiener process.
rm(list=ls()) #Removes all items in Environment!
# A typical Wiener process.
TT <- 5
N <- 10000
dt <- TT / N
Time <- seq(from = dt, to = TT, by = dt)
set.seed(560746) # for reproducibility purposes.
delta_z <- rnorm(N, 0, 1) * dt^0.5
z <- 25 + cumsum(delta_z) # The initial price is 25.Graphically:
plot(Time, z, type = "l",
ylab = "z in the texbook, it could be a stock price",
xlab = "Time, this could be 5 days or even shorter periods")
lines(Time, seq(25, (25 + sum(delta_z)), length.out = N), lty = 2,
col = "purple")
points(Time[1], 25, col = "blue", cex = 2, pch = 16)
points(Time[N], z[N], col = "red", cex = 2, pch = 16)
legend("topleft", legend = c(round(z[N], 2)), bg = "white",
text.col = "red")
Figure 4.1: Wiener process 10,000 steps from blue z(0)=25 to z(T)=26.2 in red.
This is supposed to mimic the evolution of a stock price. Now, we propose to produce not only one but 100 processes.
# 100 Wiener processes.
set.seed(365633)
mat.delta_z <- matrix(rnorm(100 * N, mean = 0, sd = 1), 100, N)
mat.z1 <- 25 + t(apply(mat.delta_z, 1, cumsum))
mat.z <- t(cbind(matrix(25, 100), mat.z1))
Time2 <- c(0, Time)
pred <- mean(mat.z[(N + 1), ])
# 95% range.
rangetop <- pred + qnorm(0.975) * var(mat.z[(N + 1), ])^0.5
rangedown <- pred - qnorm(0.975) * var(mat.z[(N + 1), ])^0.5Graphically.
matplot(Time2, mat.z, type = "l", lty = 1, col = rgb(0, 0, 1, 0.3),
ylab = "z in the texbook, it could be a stock price",
xlab = "Time")
abline(h = pred, lty = 2, col = "red")
abline(v = 0, lty = 2)
lines(Time2, seq(25, pred, length.out = N + 1), lty = 2)
points(Time2[N + 1], pred, col = "red", cex = 2, pch = 16)
points(0, 25, col = "black", cex = 2, pch = 16)
legend("topleft", legend = c(round(pred, 2)), bg = "white",
text.col = "red")
abline(h = rangetop, lty = 2)
abline(h = rangedown, lty = 2)
Figure 4.2: 100 Wiener process 10,000 steps, initial price at 25, red is the mean at T.
Let’s confirm this 95% confidence interval.
sum(mat.z[(N + 1), ] > rangetop)## [1] 3sum(mat.z[(N + 1), ] < rangedown)## [1] 2This means that 5 values are outside this 95% confidence interval. Confidence intervals are then correct.
summary(mat.z[(N + 1), ])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -155.23 -29.14 23.33 25.95 72.21 310.30A potential problem is that the values can be negative. This is problematic when we are interested to model stock prices.
Let’s explore deeper the stochastic process below.
set.seed(1)
delta_z <- rnorm(N, 0, 1) * dt^0.5
z <- 25 + cumsum(delta_z)
### A stochastic process
par(mfrow = c(2, 2), mai = c(0.4, 0.4, 0.4, 0.4))
par(pty = "s")
plot(Time, z, type = "l", xlim = c(0, TT), ylim = c(22.5, 26),
main = "z", ylab = "")
plot(Time, delta_z, type = "h", xlim = c(0, TT),
main = "delta_z", ylab = "")
abline(h = 0, lty = 2, col = "orange", lwd = 2)
plot(Time, z, type = "l", xlim = c(0, 1/12), ylim = c(24.9, 25.35),
main = "z (first month only)", ylab = "")
plot(Time, delta_z, type = "h", xlim = c(0, 1/12),
col = ifelse(delta_z < 0, "red", "blue"), ylim = c(-0.055, 0.055),
main = "delta_z (first month only)", ylab = "")
Figure 4.3: A stochastic process.
These processes are driven by a random component. The lower panel is revealing because we can see how the random component goes from positive to negative without any pattern, it is entirely random. Now, let’s explore the role of \(\Delta t\).
par(mfrow=c(2, 2), mai = c(0.4, 0.4, 0.4, 0.4))
par(pty = "s")
plot(Time[(seq(from = 1, to = N, length.out = 10))],
z[(seq(from = 1, to = N, length.out = 10))], type = "l",
ylim = c(22.5, 26), main = "Relatively large value of delta_t",
ylab = "")
abline(v = 0, lty = 2)
abline(v = 5, lty = 2)
plot(Time[(seq(from = 1, to = N, length.out = 100))],
z[(seq(from = 1, to = N, length.out = 100))], type = "l",
ylim = c(22.5, 26), ylab = "z", main = "Smaller value of delta_t")
abline(v = 0, lty = 2)
abline(v = 5 , lty = 2)
plot(Time[(seq(from = 1, to = N, length.out = 1000))],
z[(seq(from = 1, to = N, length.out = 1000))], type = "l",
ylim = c(22.5, 26), ylab = "z", main = "Smaller value of delta_t")
abline(v = 0, lty = 2)
abline(v = 5 , lty = 2)
plot(Time, z, type = "l", ylim = c(22.5, 26), # length 10,000
main = "True process as delta_t tends to 0")
abline(v = 0, lty = 2)
abline(v = 5, lty = 2)
Figure 4.4: Wiener process.
The \(\Delta t\) value determines how frequent are the changes.
A few properties of stochastic processes.
mean(delta_z) # should be zero## [1] -0.0001461726var(delta_z)^0.5 # should be dt^.5## [1] 0.02263698dt^0.5## [1] 0.02236068var(delta_z) # should be dt## [1] 0.0005124328dt## [1] 5e-044.2 The generalized Wiener process.
Now let’s analyze the case of the generalized Wiener process.
# Figure 14.2
TT <- 50
N <- 200
dt <- TT / N
Time <- seq(from = 0, to = TT, by = dt)
set.seed(123)
delta_z <- 1.5 * rnorm(N, 0, 1) * dt^0.5
delta_x <- (0.3 * dt) + delta_z
x <- c(0, cumsum(delta_x))
z <- c(0, cumsum(delta_z))Graphically:
plot(Time, x, type = "l", lwd = 2, ylim = c(-7, 22))
lines(Time, z, col = "red", lwd = 2)
lines(Time, 0.3 * Time, col = "blue", lwd = 2)
abline(0, 0)
abline(v = 0)
legend("topleft", legend = c("Generalized Wiener process",
"Drift", "Basic Wiener process"),
col = c("black", "blue", "red"), lwd = 3, bg = "white", cex = 0.8)
Figure 4.5: In Hull: Figure 14.2
Here, the generalized Wiener process is decomposed into the drift and the basic Wiener process. A zoom of the same plot.
# Figure 14.2 Zoom
plot(Time, x, type = "b", ylim = c(-0.6, 1), xlim = c(0, 1), lwd = 2)
lines(Time, z, col = "red", lwd = 2, type = "b")
lines(Time, 0.3 * Time, col = "blue", lwd = 2, type = "b")
abline(0, 0)
abline(v = 0)
abline(v = dt, lty = 2)
abline(v = dt * 2, lty = 2)
abline(v = dt * 3, lty = 2)
abline(v = dt * 4, lty = 2)
points(dt, 0, pch = 1, col = "blue", lwd = 2)
points(dt * 2, 0, pch = 1, col = "blue", lwd = 2)
points(dt * 3, 0, pch = 1, col = "blue", lwd = 2)
points(dt * 4, 0, pch = 1, col = "blue", lwd = 2)
legend("topleft", legend = c("Generalized Wiener process",
"Drift", "Basic Wiener process"),
col = c("black", "blue", "red"), lty = 1, bg = "white", lwd = 2)
Figure 4.6: In Hull: Figure 14.2 (zoom)
This is a nice visual representation of: generalized = drift + Wiener.
See how the process changes when we consider a lower value of \(a\).
set.seed(123)
delta_xlow <- (0.15 * dt) + delta_z
xlow <- c(0, cumsum(delta_xlow))
# Now plot.
plot(Time, xlow, type = "l", ylim = c(-7, 22), lwd = 2, ylab = "x")
lines(Time, x, lwd = 2, col = "grey")
lines(Time, z, lwd = 2, col = "red")
lines(Time, 0.15 * Time, lty = 2, lwd = 2)
lines(Time, 0.3 * Time, lty = 2, lwd = 2, col = "grey")
abline(0, 0)
abline(v = 0)
legend("topleft", legend = c("Original generalized Wiener process",
"Original drift", "New generalized Wiener process", "New drift",
"Basic Wiener process"), lty = c(1, 2, 1, 2, 1), bg = "white", lwd = 2,
col = c("grey", "grey", "black", "black", "red"), cex = 0.7)
Figure 4.7: Generalized Wiener process: lower a
And a higher value of \(b\).
set.seed(123)
delta_zhigh <- 3 * rnorm(N, 0, 1) * dt^0.5
delta_xhigh <- (0.3 * dt) + delta_zhigh
xhigh <- c(0, cumsum(delta_xhigh))
zhigh <- c(0, cumsum(delta_zhigh))
plot(Time, xhigh, type = "l", ylim = c(-7, 22), lwd = 2, ylab = "x")
lines(Time, x, lwd = 2, col = "grey")
lines(Time, zhigh, lwd = 2, col = "red")
lines(Time, z, lwd = 2, col = "grey")
lines(Time, 0.3 * Time, lty = 2, lwd = 2)
abline(0, 0)
abline(v = 0)
legend("topleft", legend = c("New generalized Wiener process",
"Original generalized Wiener process", "Drift" ,
"New basic Wiener process", "Original basic Wiener process"),
lty = c(1, 1, 2, 1, 1), bg = "white", lwd = 2,
col = c("black", "grey", "black", "red", "grey"), cex = 0.7)
Figure 4.8: Generalized Wiener process: higher b
Note the changes are now more pronounced.
# Table 14.1
set.seed(19256)
delta_S <- rep(0, 10)
S <- rep(100, 10)
epsilon <- rep(0, 10)
for(i in 1:10){
epsilon[i] <- rnorm(1, 0, 1)
delta_S[i] <- 0.15 * (1 / 52) * S[i] + 0.3 * ((1 / 52)^0.5) *
epsilon[i] * S[i]
S[i+1] <- S[i] + delta_S[i]
}
epsilon <- c(epsilon, NA)
delta_S <- c(delta_S, NA)
results <- data.frame(S, epsilon, delta_S)
kable(results, caption = "Simulation of stock price (Table 14.1 in Hull).",
digits = 4) |>
kable_styling(latex_options = "HOLD_position")| S | epsilon | delta_S |
|---|---|---|
| 100.0000 | 0.8269 | 3.7284 |
| 103.7284 | 0.4958 | 2.4388 |
| 106.1672 | -0.9314 | -3.8075 |
| 102.3597 | -0.1496 | -0.3418 |
| 102.0179 | 0.7310 | 3.3969 |
| 105.4148 | -0.6111 | -2.3761 |
| 103.0387 | 1.2491 | 5.6516 |
| 108.6904 | -0.6243 | -2.5096 |
| 106.1808 | 0.6995 | 3.3964 |
| 109.5772 | 0.3613 | 1.9629 |
| 111.5401 | NA | NA |
It is interesting to note that in Hull the final simulated price is 111.54 and here it is 111.5401. Is this a pure coincidence?
A correlated process.
# Section 14.5
rho <- 0.8
TT <- 1
steps <- 1000
dt <- TT / steps
Time <- seq(from = dt, to = TT, by = dt)
set.seed(123)
e1 <- rnorm(steps, 0, 1)
e2 <- rho * e1 + ((1 - rho^2)^0.5) * rnorm(steps, 0, 1)
delta_z1 <- 1.5 * e1 * dt^0.5
delta_z2 <- 1.5 * e2 * dt^0.5
delta_x1 <- (0.3 * dt) + delta_z1
delta_x2 <- (0.3 * dt) + delta_z2
x1 <- cumsum(delta_x1)
x2 <- cumsum(delta_x2)
x <- c(x1, x2)Visually:
plot(Time, x1, type = "l", ylim = c(min(x), max(x)), ylab = "x1 and x2")
lines(Time, x2, col = "red")
abline(0, 0)
Figure 4.9: Correlated processes x1 and x2.
The following code needs revision.
# Figure 14.2 with Itô process
TT <- 50
N <- 200
dt <- TT / N
Time <- seq(from = 0, to = TT, by = dt)
set.seed(1)
axt <- 0.3 + rnorm(N + 1, 0, 1) * 0.001 * Time
bxt <- 0.001 * Time + 0.1
delta_z <- bxt * rnorm(N + 1, 0, 1) * dt^0.5
delta_x <- (axt * dt) + delta_z
xI <- cumsum(delta_x)
zI <- cumsum(delta_z)
# Now plot.
plot(Time, xI, type = "l", lwd = 2, ylim = c(0, 18), ylab = "x")
lines(Time, zI, col = "red", lwd = 2)
lines(Time, axt * Time, col = "blue", lwd = 2)
abline(0, 0)
abline(v = 0)
legend("topleft", legend = c("Itô process",
"Stochastic drift", "Wiener process"),
col = c("black", "blue", "red"), lty = 1, bg = "white", lwd = 2)
Figure 4.10: Itô process. This figure needs revision.
This is a clear view of a typical Itô process. Taken from Esteban Moro.
Figure 4.11: Itô process.
The model of stock price behavior implies that a stock’s price at time \(T\), given its price today, is lognormally distributed. The following is similar to example 15.1 in Hull.
set.seed(10101010)
for(j in 1:100) {
delta_S <- NULL
S <- rep(40, 365) # Now it is 40 for all days, the loop fill this out.
epsilon <- NULL
for(i in 1:365) {
epsilon[i] <- rnorm(1, 0, 1)
delta_S[i] <- 0.16 * (1 / 365) * S[i] + 0.2 * ((1 / 365)^0.5) *
epsilon[i] * S[i]
S[i+1] <- S[i] + delta_S[i]
}
if (j==1) {
plot(S, type = "l", ylim = c(20, 80), ylab = "$", xlab = "Time (days)")
abline(h = 10, col = "red", lwd = 3)
}
lines(S)
}
abline(h = 32.55, lwd = 2, col = "red")
abline(h = 56.56, lwd = 2, col = "red")
abline(v = 180, lwd = 2, col = "red")
Figure 4.12: A geometric Brownian motion simulation: 100 stock price paths.
The code above can be problematic for a few reasons. Can you see why? An alternative code is the following. Again, this is similar to example 15.1 in Hull.
S0 <- 40
nDays <- 360
mu <- 0.16
sig <- 0.2
TT <- 0.5
mean.logST <- log(S0)+(mu - (sig^2 / 2)) * TT
variance.logST <- sig^2 * TT
lower.ST <- exp(mean.logST + qnorm(0.025) * variance.logST^0.5)
upper.ST <- exp(mean.logST + qnorm(0.975) * variance.logST^0.5)
set.seed(3)
nSim <- 1000
SP_sim <- matrix(0, nrow = nDays, ncol = nSim)
for(i in 1:nSim){
SVec <- rep(0, nDays)
SVec[1] <- S0
for(j in 2:nDays){
DeltaS <- mu * SVec[j - 1] * (1 / nDays) + sig * SVec[j - 1] *
rnorm(1) * (1 / nDays)^0.5
SVec[j] <- SVec[j - 1] + DeltaS
}
SP_sim[, i] <- SVec
}Now, let’s plot the results.
matplot(SP_sim, type = 'l', col = rgb(0, 0, 1 ,0.3), lty = 1,
ylab = 'This is similar as example 15.1',
xlab = "Time")
points(0, S0, col = "black", pch = 19)
abline(h = lower.ST)
abline(h = upper.ST)
abline(v = 180)
Figure 4.13: A geometric Brownian motion simulation: 1000 stock price paths.
Let’s verify that the there is a 95% probability that the stock price in 6 months will lie between 32.55 and 56.56.
1 - (sum(SP_sim[180, ] < lower.ST) + sum(SP_sim[180, ] > upper.ST)) / 1000## [1] 0.952S0 <- 40
mu <- 0.16
sig <- 0.2
TT <- 0.5
mean.logST <- log(S0)+(mu - (sig^2 / 2)) * TT
variance.logST <- sig^2 * TT
lower.ST <- exp(mean.logST + qnorm(0.025) * variance.logST^0.5)
upper.ST <- exp(mean.logST + qnorm(0.975) * variance.logST^0.5)
lower.ST## [1] 32.51491upper.ST## [1] 56.6029S0 <- 40
mu <- 0.16
sig <- 0.2
TT <- 1
mean.logST <- log(S0)+(mu - (sig^2 / 2)) * TT
variance.logST <- sig^2 * TT
lower.ST <- exp(mean.logST + qnorm(0.05) * variance.logST^0.5)
upper.ST <- exp(mean.logST + qnorm(0.95) * variance.logST^0.5)
lower.ST## [1] 33.11243upper.ST## [1] 63.93393Nice.
This is how a stock price evolves according to the Black-Scholes formula.
5 VaR.
Here we develop and extend the example called “Investment in Four Stock Indices”.
5.1 Prepare the data.
Load the data and prepare the database.
rm(list=ls())
df <- read.table("dataR.txt", sep = "", header = TRUE)
df.hoy <- read.table("hoy.txt", sep = "", header = TRUE)
library(zoo)
# Easy to remember names:
DJIA <- df[1]
FTSE100 <- df[2]
USDGBP <- df[3]
CAC40 <- df[4]
EURUSD <- df[5]
Nikkei <- df[6]
YENUSD <- df[7]
# Easy to remember names:
DJIA.hoy <- df.hoy[1]
FTSE100.hoy <- df.hoy[2]
USDGBP.hoy <- df.hoy[3]
CAC40.hoy <- df.hoy[4]
EURUSD.hoy <- df.hoy[5]
Nikkei.hoy <- df.hoy[6]
YENUSD.hoy <- df.hoy[7]Replicate Table 22.2 in Hull.
# Adjustments for exchange rate risk. See Hull for details.
# Now the variables starting with "A" stands for "adjusted".
# We do not need to adjust DJIA since we are supposed to be in the US.
AFTSE100 <- FTSE100 * USDGBP
ACAC40 <- CAC40 / EURUSD
ANikkei <- Nikkei / YENUSD
# Inspection of the original data.
# Verify this is the same as Table 22.2 in Hull.
Table22.2 <- data.frame(DJIA, AFTSE100, ACAC40, ANikkei)
kable(head(Table22.2), caption = "US dollar equivalent of stock indices for
historical simulation (first values).", digits = 4,
row.names = TRUE) |>
kable_styling(latex_options = "HOLD_position")| DJIA | FTSE100 | CAC40 | Nikkei | |
|---|---|---|---|---|
| 1 | 11219.38 | 11131.84 | 6373.894 | 131.7744 |
| 2 | 11173.59 | 11096.28 | 6378.162 | 134.3818 |
| 3 | 11076.18 | 11185.35 | 6474.040 | 135.9433 |
| 4 | 11124.37 | 11016.71 | 6357.486 | 135.4381 |
| 5 | 11088.02 | 11040.73 | 6364.765 | 134.1003 |
| 6 | 11097.87 | 11109.50 | 6431.668 | 136.1710 |
Now the last part of Table 22.2.
kable(tail(Table22.2), caption = "US dollar equivalent of stock indices for
historical simulation (last values).", digits = 4) |>
kable_styling(latex_options = "HOLD_position")| DJIA | FTSE100 | CAC40 | Nikkei | |
|---|---|---|---|---|
| 496 | 10609.66 | 8820.214 | 5677.136 | 112.1698 |
| 497 | 11019.69 | 8878.184 | 5689.850 | 109.5471 |
| 498 | 11388.44 | 9734.020 | 6230.006 | 111.6185 |
| 499 | 11015.69 | 9656.261 | 6181.953 | 113.2290 |
| 500 | 10825.17 | 9438.580 | 6033.935 | 114.2604 |
| 501 | 11022.06 | 9599.898 | 6200.396 | 112.8221 |
We are interested to evaluate what happened on September 26, 2008. This is why we need the real values.
AFTSE100.hoy <- FTSE100.hoy * USDGBP.hoy
ACAC40.hoy <- CAC40.hoy / EURUSD.hoy
ANikkei.hoy <- Nikkei.hoy / YENUSD.hoy
# Table.
Table22.2hoy <- data.frame(DJIA.hoy, AFTSE100.hoy, ACAC40.hoy, ANikkei.hoy)
kable(tail(Table22.2hoy), caption = "Here we show 501 and 502 values.",
digits = 4) |>
kable_styling(latex_options = "HOLD_position")| DJIA | FTSE100 | CAC40 | Nikkei |
|---|---|---|---|
| 11022.06 | 9599.898 | 6200.396 | 112.8221 |
| 11143.13 | 9379.123 | 6081.478 | 112.1679 |
# S stands for scenarios for September 26th 2008:
SDJIA <- DJIA[501, ] * DJIA[2:501, ] / DJIA[1:500, ]
SFTSE100 <- AFTSE100[501, ] * AFTSE100[2:501, ] / AFTSE100[1:500, ]
SCAC40 <- ACAC40[501, ] * ACAC40[2:501, ] / ACAC40[1:500, ]
SNikkei <- ANikkei[501, ] * ANikkei[2:501, ] / ANikkei[1:500, ]
# Real values for September 26th 2008:
SDJIA.hoy <- DJIA.hoy[2, ] * DJIA.hoy[2, ] / DJIA.hoy[1, ]
SFTSE100.hoy <- AFTSE100.hoy[2, ] * AFTSE100.hoy[2, ] / AFTSE100.hoy[1, ]
SCAC40.hoy <- ACAC40.hoy[2, ] * ACAC40.hoy[2, ] / ACAC40.hoy[1, ]
SNikkei.hoy <- ANikkei.hoy[2, ] * ANikkei.hoy[2, ] / ANikkei.hoy[1, ]
# Investment portfolio value
p <- (4000 * DJIA[2:501, ] / DJIA[1:500, ]) +
(3000 * AFTSE100[2:501, ] / AFTSE100[1:500, ]) +
(1000 * ACAC40[2:501, ] / ACAC40[1:500, ]) +
(2000 * ANikkei[2:501,] / ANikkei[1:500, ])
# Evaluated according to the real values.
p.hoy <- (4000 * DJIA.hoy[2, ] / DJIA.hoy[1, ]) +
(3000 * AFTSE100.hoy[2, ] / AFTSE100.hoy[1, ]) +
(1000 * ACAC40.hoy[2, ] / ACAC40.hoy[1, ]) +
(2000 * ANikkei.hoy[2, ] / ANikkei.hoy[1, ])
# Losses l
l <- 10000 - p
l.hoy <- 10000 - p.hoyReplicate Table 21.3 in Hull.
# Verify this is the same as Table 21.3 in Hull.
Table21.3 <- data.frame(SDJIA, SFTSE100, SCAC40, SNikkei, p, l)
kable(head(Table21.3), caption = "Scenarios generated for September 26,
2008 (first values).", digits = 4, row.names = TRUE) |>
kable_styling(latex_options = "HOLD_position")| SDJIA | SFTSE100 | SCAC40 | SNikkei | p | l | |
|---|---|---|---|---|---|---|
| 1 | 10977.08 | 9569.230 | 6204.547 | 115.0545 | 10014.334 | -14.3338 |
| 2 | 10925.97 | 9676.957 | 6293.602 | 114.1331 | 10027.481 | -27.4813 |
| 3 | 11070.01 | 9455.160 | 6088.768 | 112.4028 | 9946.736 | 53.2641 |
| 4 | 10986.04 | 9620.831 | 6207.495 | 111.7077 | 9974.861 | 25.1394 |
| 5 | 11031.85 | 9659.698 | 6265.572 | 114.5642 | 10063.635 | -63.6351 |
| 6 | 11153.55 | 9660.357 | 6305.745 | 112.9480 | 10085.834 | -85.8336 |
The last values of Table 21.3.
kable(tail(Table21.3), caption = "Scenarios generated for September 26,
2008 (last values).", digits = 4) |>
kable_styling(latex_options = "HOLD_position")| SDJIA | SFTSE100 | SCAC40 | SNikkei | p | l | |
|---|---|---|---|---|---|---|
| 495 | 10574.20 | 9495.783 | 6099.845 | 114.3359 | 9815.550 | 184.4496 |
| 496 | 11448.03 | 9662.993 | 6214.283 | 110.1841 | 10129.781 | -129.7806 |
| 497 | 11390.89 | 10525.305 | 6789.019 | 114.9555 | 10555.795 | -555.7954 |
| 498 | 10661.30 | 9523.211 | 6152.571 | 114.4499 | 9866.256 | 133.7444 |
| 499 | 10831.43 | 9383.488 | 6051.936 | 113.8498 | 9857.465 | 142.5355 |
| 500 | 11222.53 | 9763.974 | 6371.450 | 111.4019 | 10126.439 | -126.4390 |
Nice. Values are exactly the same as in Hull.
See what really happened.
Table21.3hoy <- data.frame(SDJIA.hoy, SFTSE100.hoy, SCAC40.hoy, SNikkei.hoy,
p.hoy,l.hoy)
kable(head(Table21.3hoy), caption = "What really happened.", digits = 4,
row.names = TRUE) |>
kable_styling(latex_options = "HOLD_position")| SDJIA.hoy | SFTSE100.hoy | SCAC40.hoy | SNikkei.hoy | p.hoy | l.hoy | |
|---|---|---|---|---|---|---|
| 1 | 11265.53 | 9163.425 | 5964.841 | 111.5174 | 9944.167 | 55.8325 |
Visually.
plot(l, type = "h", lwd = 2,
xlab = "500 scenarios (sorted by historical date)",
ylab = "Gains (-) and losses (+)", col = ifelse(l < 0, "blue", "red"))
points(501, l.hoy, pch = 19, cex = 2)
abline(0, 0)
abline(v = 0)
lines(seq(0, max(l), length.out = 500), lty = 2)
lines(seq(0, min(l), length.out = 500), lty = 2)
abline(h = 253.385, lwd = 2, col = "green")
legend("bottomleft", legend = c("Gains", "Losses", "Historic VaR (253.385)"),
col = c("blue", "red", "green"), lty = 1, bg = "white", lwd = 2)
Figure 5.1: Today is September 25, 2008. What is the worst that can happen to my portfolio tomorrow?
Now, a histogram of losses.
hist(l, 15, xlab = "Losses (+) and gains (-)", main = NULL)
points(l.hoy, 1, pch = 19, cex = 2)
abline(v = 253.385, col = "red", lwd = 2)
legend("topleft", legend = c("Historic VaR (253.385)"), col = "red",
bg = "white", lwd = 2)
Figure 5.2: Histogram of losses for the scenarios considered between September 25 and 26, 2008: a histogram.
The point refers to what really happened, and the red line is the VaR following a historic approach. Now, the same information but now using a density plot.
densi <- density(l)
plot(densi, main = "", xlab = "Losses (+) and gains (-)")
polygon(densi, col = "grey", border = "black")
points(l.hoy, 0, pch = 19, cex = 2)
abline(v = 253.385, col = "red", lwd = 2)
legend("topleft", legend = c("Historic VaR (253.385)"), col = "red",
bg = "white", lwd = 2)
Figure 5.3: Histogram of losses for the scenarios considered between September 25 and 26, 2008: a density plot.
Looks nice.
Here is a histogram that looks like the one in Hull.
hist(l, 15, axes = F, lty = "blank", col = "grey", main = "",
xlab = "Losses (+) and gains (-)")
axis(1, pos = 0)
axis(2, pos = 0, las = 1)
Figure 5.4: Figure 21.3. Histogram of losses for the scenarios considered between September 25 and 26, 2008.
A table for the losses ranked.
Table22.4 <- data.frame(scenario.number = order(l, decreasing = TRUE),
`loss in positive values` =
l[order(l, decreasing = TRUE)])
kable(head(Table22.4), caption = "Losses ranked from highest to lowest for
500 scenarios (first values).", digits = 4, row.names = TRUE) |>
kable_styling(latex_options = "HOLD_position")| scenario.number | loss.in.positive.values | |
|---|---|---|
| 1 | 494 | 477.8410 |
| 2 | 339 | 345.4351 |
| 3 | 349 | 282.2038 |
| 4 | 329 | 277.0413 |
| 5 | 487 | 253.3850 |
| 6 | 227 | 217.9740 |
And the last values (not shown in Hull).
kable(tail(Table22.4), caption = "Losses ranked from highest to lowest for
500 scenarios (last values).", digits = 4) |>
kable_styling(latex_options = "HOLD_position")| scenario.number | loss.in.positive.values | |
|---|---|---|
| 495 | 395 | -224.5146 |
| 496 | 489 | -284.9247 |
| 497 | 341 | -307.9301 |
| 498 | 377 | -316.4893 |
| 499 | 379 | -333.0218 |
| 500 | 497 | -555.7954 |
Note that one of the highest gains happen in scenario 497, almost at the end.
# Results
(VaR.hist <- Table22.4[5, 2])## [1] 253.385The VaR according to the historic approach is the highest 5th loss because the top 1% of 500 scenarios is the 5th scenario. We can view this in a barplot.
barplot(l[order(l, decreasing = TRUE)],
names.arg = order(l, decreasing = TRUE), ylim = c(0, 500),
xlim = c(0, 20), las = 2, cex.names = 1, ylab = "Loss ($000s)",
col = "red", xlab = "Scenario number (ranked)")
abline(h = 253.385, lwd = 3)
abline(h = l.hoy, lwd = 3, col = "green")
legend("topright", legend = c("Historic VaR (253.385)",
"What really happened"), col = c("black", "green"),
bg = "white", lwd = 2)
Figure 5.5: Extract of table 22.4 in a graph VaR historical approach 1 percent 253.385 What really happen in green.
quantile(l, 0.99)## 99%
## 218.3281barplot(quantile(l), ylim = c(min(l), max(l)),
col = ifelse(quantile(l) < 0, "blue", "red"),
xlab = "Quantiles", ylab = "Gains (-) and losses (+)")
Figure 5.6: Ilustration of portfolio losses in quantiles.
# quantile(l,seq(from=0.8,to=0.99,by=0.001))# Results
VaR.hist.quantile <- quantile(l, 0.99)
kable(VaR.hist.quantile, caption = "Quantile.", digits = 4) |>
kable_styling(latex_options = "HOLD_position")| x | |
|---|---|
| 99% | 218.3281 |
plot(seq(from = 0.902, to = 1, by = 0.002),
quantile(l, seq(from = 0.902, to = 1, by = 0.002)), cex = 1.5,
xlab = "Confidence level (vertical line at 99%)", ylab = "VaR")
abline(v = 0.99, lty = 2, lwd = 2)
abline(h = quantile(l, 0.99), col = "blue", lty = 2, lwd = 2)
abline(h = Table22.4[5, 2], col = "red", lty = 2, lwd = 2)
legend("topleft", legend = c("99% level", "Historic VaR (253.385)",
"99% quantile VaR (218.3281)"),
col = c("black", "red", "blue"), lty = 2, bg = "white", lwd = 2)
Figure 5.7: If we use a quantile approach.
Summary of results.
Results <- data.frame(VaR.method = c("VaR.hist","VaR.hist.quantile"),
value = c(VaR.hist,VaR.hist.quantile))
kable(Results, caption = "Results.", digits = 4) |>
kable_styling(latex_options = "HOLD_position")| VaR.method | value | |
|---|---|---|
| VaR.hist | 253.3850 | |
| 99% | VaR.hist.quantile | 218.3281 |
5.2 Model building approach.
Let’s calculate more VaR with different methods.
# Model building example -- equally weighted
# Compute returns (this code should be simplified)
RDJIA <- (DJIA[2:501, ] - DJIA[1:500, ]) / DJIA[1:500, ]
RAFTSE100 <- (AFTSE100[2:501, ] - AFTSE100[1:500, ]) / AFTSE100[1:500, ]
RACAC40 <- (ACAC40[2:501, ] - ACAC40[1:500, ]) / ACAC40[1:500, ]
RANikkei <- (ANikkei[2:501, ] - ANikkei[1:500, ]) / ANikkei[1:500, ]
# Gather returns.
R <- merge.zoo(RDJIA, RAFTSE100, RACAC40, RANikkei)
# Compute returns for "today".
RDJIA.hoy <- (DJIA.hoy[2, ] - DJIA.hoy[1, ]) / DJIA.hoy[1, ]
RAFTSE100.hoy <- (AFTSE100.hoy[2, ] - AFTSE100.hoy[1, ]) / AFTSE100.hoy[1, ]
RACAC40.hoy <- (ACAC40.hoy[2, ] - ACAC40.hoy[1, ]) / ACAC40.hoy[1, ]
RANikkei.hoy <- (ANikkei.hoy[2, ] - ANikkei.hoy[1, ]) / ANikkei.hoy[1, ]
# Gather returns.
R.hoy <- merge.zoo(RDJIA.hoy, RAFTSE100.hoy, RACAC40.hoy, RANikkei.hoy)Replicate Hull’s table 23.5.
# correlation
Rcor <- cor(R)
Table23.5 <- Rcor
kable(Table23.5, caption = "Correlation matrix on September 25, 2008,
calculated by giving equal weight to the last 500 daily returns.",
digits = 4) |>
kable_styling(latex_options = "HOLD_position")| RDJIA | RAFTSE100 | RACAC40 | RANikkei | |
|---|---|---|---|---|
| RDJIA | 1.0000 | 0.4891 | 0.4957 | -0.0619 |
| RAFTSE100 | 0.4891 | 1.0000 | 0.9181 | 0.2009 |
| RACAC40 | 0.4957 | 0.9181 | 1.0000 | 0.2110 |
| RANikkei | -0.0619 | 0.2009 | 0.2110 | 1.0000 |
Replicate Hull’s table 23.6.
# covariance
Rcov <- cov(R)
Table23.6 <- Rcov
kable(Table23.6, caption = "Covariance matrix on September 25, 2008,
calculated by giving equal weight to the last 500 daily returns.",
digits = 8) |>
kable_styling(latex_options = "HOLD_position")| RDJIA | RAFTSE100 | RACAC40 | RANikkei | |
|---|---|---|---|---|
| RDJIA | 0.00012295 | 0.00007697 | 0.00007683 | -0.00000949 |
| RAFTSE100 | 0.00007697 | 0.00020140 | 0.00018211 | 0.00003944 |
| RACAC40 | 0.00007683 | 0.00018211 | 0.00019535 | 0.00004078 |
| RANikkei | -0.00000949 | 0.00003944 | 0.00004078 | 0.00019131 |
This covariance matrix is quite similar to the one in Hull. We are not taking round values this explains minor differences.
This covariance matrix gives the variance of the portfolio losses ($000s) as
# standard deviation
Rsd <- apply(R, 2, sd)
# allocation of the 10,000 usd
alpha <- c(4000, 3000, 1000, 2000)
# Answers can be compared in page 538 (Chapter 23)
(Pvar <- t(alpha) %*% Rcov %*% alpha) # equation in page 506 (chapter 22)## [,1]
## [1,] 8779.392This is slightly different as in Hull: 8,761.833. The standard deviation is the square root of this
(Psd <- Pvar^0.5)## [,1]
## [1,] 93.69841As expected, this is a bit different as in Hull (93.60). The one-day 99% VaR in is therefore
(VaR.ew <- qnorm(0.99, 0, 1) * Psd) # I did not use the simplification 2.33## [,1]
## [1,] 217.9751In Hull, this value is 217,757, which compares with 253,385, calculated using the historical simulation approach in Section 22.2.
Results <- data.frame(VaR.method = c("VaR.hist", "VaR.hist.quantile",
"VaR.ew"),
value = c(VaR.hist, VaR.hist.quantile, VaR.ew))
kable(Results, caption = "Summary of results.", digits = 4) |>
kable_styling(latex_options = "HOLD_position")| VaR.method | value |
|---|---|
| VaR.hist | 253.3850 |
| VaR.hist.quantile | 218.3281 |
| VaR.ew | 217.9751 |
We are interested to compare the equally weighted approach with the EWMA (exponentially weighted moving average method). In fact, there are two EWMA versions here, one was taken from a different author, and the other is mine (Martín).
First the other author’s approach to calculate an EWMA covariance considering \(\lambda=0.94\).
# Model building example -- EWMA
covEWMA <-
function(factors, lambda = 0.96, return.cor = FALSE) {
factor.names = colnames(factors)
t.factor = nrow(factors)
k.factor = ncol(factors)
factors = as.matrix(factors)
t.names = rownames(factors)
cov.f.ewma = array(,c(t.factor, k.factor, k.factor))
cov.f = var(factors) # unconditional variance as EWMA at time = 0
FF = (factors[1, ] - mean(factors)) %*% t(factors[1,] - mean(factors))
cov.f.ewma[1,,] = (1 - lambda) * FF + lambda * cov.f
for (i in 2:t.factor) {
FF = (factors[i, ] - mean(factors)) %*% t(factors[i, ] - mean(factors))
cov.f.ewma[i,,] = (1 - lambda) * FF + lambda * cov.f.ewma[(i-1),,]
}
dimnames(cov.f.ewma) = list(t.names, factor.names, factor.names)
if(return.cor) {
cor.f.ewma = cov.f.ewma
for (i in 1:dim(cor.f.ewma)[1]) {
cor.f.ewma[i, , ] = cov2cor(cov.f.ewma[i, ,])
}
return(cor.f.ewma)
} else {
return(cov.f.ewma)
}
}Let’s implement this approach to our problem.
Rdf <- as.data.frame(R)
C <- covEWMA(Rdf, lambda = 0.94)[500,,] # Table 23.7
Table23.7 <- C
kable(Table23.7, caption = "Covariance matrix on September 25, 2008,
calculated using the EWMA (other author's approach).", digits = 8) |>
kable_styling(latex_options = "HOLD_position")| RDJIA | RAFTSE100 | RACAC40 | RANikkei | |
|---|---|---|---|---|
| RDJIA | 0.00047986 | 0.00043008 | 0.00042553 | -0.00003989 |
| RAFTSE100 | 0.00043008 | 0.00103126 | 0.00096296 | 0.00020934 |
| RACAC40 | 0.00042553 | 0.00096296 | 0.00095350 | 0.00016796 |
| RANikkei | -0.00003989 | 0.00020934 | 0.00016796 | 0.00025383 |
The variance of portfolio losses is:
(Pvar_EWMA <- t(alpha) %*% C %*% alpha)## [,1]
## [1,] 40977.52In Hull, this value is 40,995.765. The standard deviation is the square root of this, or:
(Psd_EWMA <- Pvar_EWMA^0.5)## [,1]
## [1,] 202.4291In Hull, this value is 202.474. The one-day 99% VaR is therefore:
(VaR.ewma <- qnorm(0.99, 0, 1) * Psd_EWMA)## [,1]
## [1,] 470.9204In Hull, this value is 471,025, over twice as high as the value given when returns are equally weighted.
Results <- data.frame(VaR.method=c("VaR.hist", "VaR.hist.quantile",
"VaR.ew","VaR.ewma"),
value = c(VaR.hist, VaR.hist.quantile, VaR.ew,VaR.ewma))
kable(Results, caption = "Summary of results.", digits = 4) |>
kable_styling(latex_options = "HOLD_position")| VaR.method | value |
|---|---|
| VaR.hist | 253.3850 |
| VaR.hist.quantile | 218.3281 |
| VaR.ew | 217.9751 |
| VaR.ewma | 470.9204 |
Let’s replicate Table 23.8.
#Table 23.8
Table23.8 <- data.frame(Rsd * 100, diag(C)^0.5 * 100)
colnames(Table23.8) <- c("Equal weighting", "EWMA")
rownames(Table23.8) <- c("DJIA", "FTSE 100", "CAC 40", "Nikkei 225")
kable(t(Table23.8), caption = "Volatilities % per day using equal
weighting and EWMA (other author approach).", digits = 4) |>
kable_styling(latex_options = "HOLD_position")| DJIA | FTSE 100 | CAC 40 | Nikkei 225 | |
|---|---|---|---|---|
| Equal weighting | 1.1088 | 1.4191 | 1.3977 | 1.3832 |
| EWMA | 2.1906 | 3.2113 | 3.0879 | 1.5932 |
The estimated daily standard deviations are much higher when EWMA is used than when data are equally weighted. This is because volatilities were much higher during the period immediately preceding September 25, 2008, than during the rest of the 500 days covered by the data.
Now EWMA with my own approach.
# Model building example -- EWMA -- INEFFICIENT BUT ACCURATE.
# My approach based on Hull.
lambda = 0.94
Var.R.EWMA <- matrix(0, 501, 4)
for(i in 1:501) {
if(i==1) {
Var.R.EWMA[, (1:4)] <- Rsd^2
} else if (i > 1) {
Var.R.EWMA[i,] <- (Var.R.EWMA[(i - 1), (1:4)] * lambda) +
(R[(i - 1), (1:4)]^2 * (1 - lambda))}
}
Cov.DF.EWMA <- NULL
Cov.DC.EWMA <- NULL
Cov.DN.EWMA <- NULL
Cov.FC.EWMA <- NULL
Cov.FN.EWMA <- NULL
Cov.CN.EWMA <- NULL
for(i in 1:501){
if(i==1){
Cov.DF.EWMA[i] <- cov(RDJIA, RAFTSE100)
Cov.DC.EWMA[i] <- cov(RDJIA, RACAC40)
Cov.DN.EWMA[i] <- cov(RDJIA, RANikkei)
Cov.FC.EWMA[i] <- cov(RAFTSE100, RACAC40)
Cov.FN.EWMA[i] <- cov(RAFTSE100, RANikkei)
Cov.CN.EWMA[i] <- cov(RACAC40, RANikkei)
} else if (i > 1) {
Cov.DF.EWMA[i] <- (Cov.DF.EWMA[i - 1] * lambda) +
(RAFTSE100[i - 1] * RDJIA[i - 1] * (1 - lambda))
Cov.DC.EWMA[i] <- (Cov.DC.EWMA[i - 1] * lambda) +
(RACAC40[i - 1] * RDJIA[i - 1] * (1 - lambda))
Cov.DN.EWMA[i] <- (Cov.DN.EWMA[i - 1] * lambda) +
(RANikkei[i - 1] * RDJIA[i - 1] * (1 - lambda))
Cov.FC.EWMA[i] <- (Cov.FC.EWMA[i - 1] * lambda) +
(RACAC40[i - 1] * RAFTSE100[i - 1] * (1 - lambda))
Cov.FN.EWMA[i] <- (Cov.FN.EWMA[i - 1] * lambda) +
(RANikkei[i-1] * RAFTSE100[i-1] * (1 - lambda))
Cov.CN.EWMA[i] <- (Cov.CN.EWMA[i - 1] * lambda) +
(RANikkei[i - 1] * RACAC40[i - 1] * (1 - lambda))
}
}
C11 <- c(Var.R.EWMA[501, 1])
C22 <- c(Var.R.EWMA[501, 2])
C33 <- c(Var.R.EWMA[501, 3])
C44 <- c(Var.R.EWMA[501, 4])
C12 <- c(Cov.DF.EWMA[501])
C13 <- c(Cov.DC.EWMA[501])
C14 <- c(Cov.DN.EWMA[501])
C23 <- c(Cov.FC.EWMA[501])
C24 <- c(Cov.FN.EWMA[501])
C34 <- c(Cov.CN.EWMA[501])
VCV <- c(C11 ,C12, C13, C14, C12, C22, C23, C24,C13, C23,
C33, C34, C14, C24, C34, C44)
VCV <- matrix(VCV, 4, 4)The variance of portfolio losses is now:
(Pvar_EWMA1 <- t(alpha) %*% VCV %*% alpha)## [,1]
## [1,] 40995.76Which is basically the same as in Hull: 40,995.765. The standard deviation is the square root of this, or:
(Psd_EWMA1 <- Pvar_EWMA1^0.5)## [,1]
## [1,] 202.4741The same as in Hull: 202.474. The one-day 99% VaR is therefore:
# My approach is longer but more accurate
(VaR.ewma.martin <- qnorm(0.99, 0, 1) * Psd_EWMA1) ## [,1]
## [1,] 471.0252The same as in Hull: 471.0252. Let’s report all results so far.
Results <- data.frame(VaR.method=
c("VaR.hist",
"VaR.hist.quantile",
"VaR.ew",
"VaR.ewma",
"VaR.ewma.martin"),
value=c(VaR.hist,
VaR.hist.quantile,
VaR.ew,
VaR.ewma,
VaR.ewma.martin))
kable(Results, caption = "Summary of results.", digits = 4) |>
kable_styling(latex_options = "HOLD_position")| VaR.method | value |
|---|---|
| VaR.hist | 253.3850 |
| VaR.hist.quantile | 218.3281 |
| VaR.ew | 217.9751 |
| VaR.ewma | 470.9204 |
| VaR.ewma.martin | 471.0252 |
Let’s replicate an extended Table 23.8.
#Table 23.8
Table23.8 <- data.frame(Rsd * 100, diag(C)^0.5 * 100, diag(VCV)^0.5 * 100)
colnames(Table23.8) <- c("Equal weighting", "EWMA", "EWMA Martín")
rownames(Table23.8) <- c("DJIA", "FTSE 100", "CAC 40", "Nikkei 225")
kable(t(Table23.8), caption = "Volatilities % per day using equal
weighting and EWMA (other author approach).", digits = 4) |>
kable_styling(latex_options = "HOLD_position")| DJIA | FTSE 100 | CAC 40 | Nikkei 225 | |
|---|---|---|---|---|
| Equal weighting | 1.1088 | 1.4191 | 1.3977 | 1.3832 |
| EWMA | 2.1906 | 3.2113 | 3.0879 | 1.5932 |
| EWMA Martín | 2.1911 | 3.2115 | 3.0880 | 1.5941 |
Again: The estimated daily standard deviations are much higher when EWMA is used than when data are equally weighted. This is because volatilities were much higher during the period immediately preceding September 25, 2008, than during the rest of the 500 days covered by the data.
5.3 Optimal allocation of resources and VaR.
The current portfolio weights are:
W_hull <- alpha / 10000
W_hullt <- data.frame(W_hull)
rownames(W_hullt) <- c("DJIA", "FTSE 100", "CAC 40", "Nikkei 225")
kable(t(W_hullt), caption = "Original (and presumible sub-optimal) Hull's
portfolio weights.", digits = 4) |>
kable_styling(latex_options = "HOLD_position")| DJIA | FTSE 100 | CAC 40 | Nikkei 225 | |
|---|---|---|---|---|
| W_hull | 0.4 | 0.3 | 0.1 | 0.2 |
Let’s calculate optimal weights.
# Further extension: propose an optimal allocation of resources
R <- R * 100
R.hoy <- R.hoy * 100
stocks<- c("DJIA", "FTSE100", "CAC40", "Nikkei")
sigma <- var(R)
sd <- diag(sigma)^0.5
E <- colMeans(R)
E.hoy <- colMeans(R.hoy)
ones <- abs(E / E)
a <- c(t(ones) %*% solve(sigma) %*% ones)
b <- c(t(ones) %*% solve(sigma) %*% E)
c <- c(t(E) %*% solve(sigma) %*% E)
d <- c(a * c - (b^2))
g <- c(solve(sigma) %*% (c * ones - b * E) / d)
h <- c(solve(sigma) %*% (a * E - b * ones) / d)
ER <- seq(from = -0.3, to = 0.3, by = 0.0001)
S <- ER
W <- matrix(0, nrow = length(ER), ncol = length(stocks))
for(i in 1:length(ER)){
W[i,] <- g + h * ER[i]
ER[i] <- W[i,] %*% E
S[i] <- (t(W[i,]) %*% sigma %*% W[i, ])^0.5
}
W_min <- (solve(sigma) %*% ones) / a
R_min <- c(t(W_min) %*% E)
R_min.realized <- c(t(W_min) %*% E.hoy)
S_min <- c(t(W_min) %*% sigma %*% W_min)^0.5
# Allocation to mimic the return of the CAC40 but with lower risk
W.cac40 <- g + h * E[3]
R_cac40 <- c(t(W.cac40) %*% E)
R_cac40.realized <- c(t(W.cac40) %*% E.hoy)
S_cac40 <- c(t(W.cac40) %*% sigma %*% W.cac40)^0.5
R_hull.realized = c(t(W_hull) %*% E.hoy)Now we have four different portfolio weights alternatives. Two of them are optimal. Can we reduce the VaR by using optimal weights?
W <- data.frame(W_hull, W_min, W.cac40, c(0.25,0.25,0.25,0.25))
rownames(W) <- c("DJIA", "FTSE 100", "CAC 40", "Nikkei 225")
colnames(W) <- c("Hull (original weights)", "minimum variance (Martín)", "CAC 40 target return", "naive (1/N)")
kable(t(W), caption = "Four different weights alternatives.",
digits = 4) |>
kable_styling(latex_options = "HOLD_position")| DJIA | FTSE 100 | CAC 40 | Nikkei 225 | |
|---|---|---|---|---|
| Hull (original weights) | 0.4000 | 0.3000 | 0.1000 | 0.2000 |
| minimum variance (Martín) | 0.5592 | 0.0379 | 0.0228 | 0.3800 |
| CAC 40 target return | 0.6075 | -0.4034 | 0.4598 | 0.3361 |
| naive (1/N) | 0.2500 | 0.2500 | 0.2500 | 0.2500 |
See the results graphically.
W_hull <- alpha / 10000
R_hull <- c(t(W_hull) %*% E)
S_hull <- c(t(W_hull) %*% sigma %*% W_hull)^0.5
plot(S, ER, type = "l", lwd = 2, xlim = c(0.7, 1.6), ylim = c(-0.022, 0.005),
ylab = "Expected return", xlab = "Standard deviation")
points(diag(sigma)^0.5, E)
abline(0, R_hull / S_hull, lty = 3, col = "blue")
abline(0, R_min / S_min, lty = 3, col = "red")
points(S_min, R_min, col = "red", pch = 19, cex = 2)
abline(0, 0, lty = 3)
text(S_min, R_min, "Martin's
alpha", pos = 2, cex = 0.8)
text(sd, E, stocks, adj = -0.5, cex = 0.8)
points(S_hull, R_hull, col = "blue", pch = 19, cex = 2)
text(S_hull, R_hull, "Hull's alpha", adj = -0.3, cex = 0.8, pos = 1)
legend("bottomleft", legend = c("DJ=0.4, FTSE=0.3, CAC=0.1, Nikkei=0.2",
"DJ=0.56, FTSE=0.038, CAC=0.022, Nikkei=0.38"),
col = c("blue","red"), pch = 19, bg = "white", cex = 0.8)
Figure 5.8: Hull’s allocation is sub-optimal, the red one is a better alternative, although with negative expected return.
The figure suggest that the minimum variance alternative reduce the original portfolio risk. In fact, the minimum variance portfolio has a higher return and lower risk. This representation is an estimate of the future. We can evaluate what really happened by using the real returns.
(R_min.realized <- c(t(W_min) %*% E.hoy))## [1] 0.2629549(R_hull.realized = c(t(W_hull) %*% E.hoy))## [1] -0.5583249W_hull <- alpha / 10000
R_hull <- c(t(W_hull) %*% E)
S_hull <- c(t(W_hull) %*% sigma %*% W_hull)^0.5
W_1N <- c(0.25,0.25,0.25,0.25)
R_1N <- c(t(W_1N) %*% E)
S_1N <- c(t(W_1N) %*% sigma %*% W_1N)^0.5
plot(S, ER, type = "l", lwd = 2, xlim = c(0.7, 1.6), ylim = c(-0.022, 0.005),
ylab = "Expected return", xlab = "Standard deviation")
points(diag(sigma)^0.5, E)
abline(0, R_hull / S_hull, lty = 3, col = "blue")
abline(0, R_min / S_min, lty = 3, col = "red")
abline(0, R_1N / S_1N, lty = 3, col = "purple")
points(S_min, R_min, col = "red", pch = 19, cex = 2)
points(S_1N, R_1N, col = "purple", pch = 19, cex = 2)
abline(0, 0, lty = 3)
text(S_min, R_min, "Martin's
alpha", pos = 2, cex = 0.8)
text(sd, E, stocks, adj = -0.5, cex = 0.8)
points(S_hull, R_hull, col = "blue", pch = 19, cex = 2)
text(S_hull, R_hull, "Hull's alpha", adj = -0.3, cex = 0.8, pos = 1)
text(S_1N, R_1N, "1/N alpha", adj = -0.3, cex = 0.8, pos = 3)
legend("bottomleft",
legend = c("Hull: DJ=0.4, FTSE=0.3, CAC=0.1, Nikkei=0.2",
"Martin: DJ=0.56, FTSE=0.038, CAC=0.022, Nikkei=0.38",
"1/N: DJ=0.25, FTSE=0.25, CAC=0.25, Nikkei=0.25"),
col = c("blue","red", "purple"), pch = 19, bg = "white", cex = 0.8)
Figure 5.9: Hull’s allocation is sub-optimal, the red one is a better alternative, although with negative expected return.
Hull’s sub-optimal portfolio allocation turns out to be problematic. His recommendation would lead to a loss and the optimal minimum variance would lead to a positive return. Nice.
set.seed <- 1
M <- rnorm(10000, R_min, S_min)
hist(M, 100, xlim = c(-4, 4), ylim = c(0, 450), main = "")
abline(v = R_min.realized, col = "blue", lwd = 3)
abline(v = quantile(M, 0.025), lty = 2)
abline(v = quantile(M, 0.975), lty = 2)
Figure 5.10: Martin recommendation and realized outcome in blue.
set.seed = 1
H = rnorm(10000, R_hull, S_hull)
hist(H,100, xlim = c(-4, 4), ylim = c(0, 450), main = "")
abline(v = R_hull.realized, col = "blue", lwd = 3)
abline(v = quantile(H, 0.025), lty = 2)
abline(v = quantile(H, 0.975), lty = 2)
Figure 5.11: Hull recommendation and realized outcome.
set.seed <- 1
Ca <- rnorm(10000, R_cac40, S_cac40)
hist(Ca, 100, xlim = c(-4, 4), ylim = c(0, 450), main = "")
abline(v = R_cac40.realized, col = "blue", lwd = 3)
abline(v = quantile(C, 0.025), lty = 2)
abline(v = quantile(C, 0.975), lty = 2)
Figure 5.12: CAC40 recommendation and realized outcome.
Let’s compare a different portfolio. This is optimal and targets the CAC 40 expected return.
plot(S, ER, type = "l", lwd = 2, xlim = c(0.7, 1.6),
ylim = c(-0.022, 0.007), ylab = "Expected return",
xlab = "Standard deviation")
points(diag(sigma)^0.5, E)
abline(0, 0, lty = 3)
points(S_cac40, R_cac40, col = "red", pch = 19, cex = 2)
points(0.68, R_cac40.realized, col = "red", pch = 15, cex = 2)
text(S_cac40, R_cac40, "CAC40
expected return", pos = 2, cex = 0.8)
text(sd,E,stocks,adj=-0.5,cex=0.8)
W_hull = alpha / 10000
R_hull = c(t(W_hull) %*% E)
S_hull = c(t(W_hull) %*% sigma %*% W_hull)^0.5
points(S_hull, R_hull, col = "blue", pch = 19, cex = 2)
points(0.68, R_hull.realized, col = "blue", pch = 15, cex = 2)
abline(0,R_hull / S_hull, lty = 3, col = "blue")
abline(0,R_cac40 / S_cac40, lty = 3, col = "red")
abline(0,E[3] / sd[3], lty = 3)
abline(v = S_cac40, lty = 3)
text(S_hull, R_hull, "Hull's alpha", adj = -0.3, cex = 0.8, pos = 1)
legend("bottomleft", legend = c("DJ=0.4, FTSE=0.3, CAC=0.1, Nikkei=0.2",
"DJ=0.607, FTSE=-0.403, CAC=0.459, Nikkei=0.336"),
col = c("blue", "red"), pch = 19, bg = "white", cex = 0.8)
Figure 5.13: Estimation of a new optimal alpha: Replicate CAC40 ER.
(R_min.realized)## [1] 0.2629549(R_hull.realized)## [1] -0.5583249(R_cac40.realized)## [1] 0.5183212As stated before, Hull’s sub-optimal portfolio allocation turns out to be problematic. His recommendation would lead to a loss (-0.558) and the optimal minimum variance would lead to a positive return (0.2629), and the optimal portfolio that targets the CAC 40 expected return leads to an even higher return (0.518). Nice.
5.4 Model building approach with optimal weights.
Now that we have two additional optimal portfolios, we can compute the corresponding VaR for the model building approach. First, the equally weighted case.
# Model building approach, equally weighted.
alpha2 <- c(W_min * 10000)
alpha3 <- c(W.cac40 * 10000)
# This is the minimum variance optimal portfolio.
Pvar2 <- t(alpha2) %*% Rcov %*% alpha2
Psd2 <- Pvar2^0.5
VaR.ew.optimal <- qnorm(0.99, 0, 1) * Psd2
# This is optimal portfolio targeting the CAC 40 expected return.
Pvar3 <- t(alpha3) %*% Rcov %*% alpha3
Psd3 <- Pvar3^0.5
VaR.ew.optimal.cac <- qnorm(0.99, 0, 1) * Psd3Summary of results:
Results <- data.frame(VaR.method=c("VaR.hist",
"VaR.hist.quantile",
"VaR.ew",
"VaR.ew.optimal",
"VaR.ew.optimal.cac",
"VaR.ewma",
"VaR.ewma.martin"),
value=c(VaR.hist,
VaR.hist.quantile,
VaR.ew,
VaR.ew.optimal,
VaR.ew.optimal.cac,
VaR.ewma,
VaR.ewma.martin))
kable(Results, caption = "Summary of results.", digits = 4) |>
kable_styling(latex_options = "HOLD_position")| VaR.method | value |
|---|---|
| VaR.hist | 253.3850 |
| VaR.hist.quantile | 218.3281 |
| VaR.ew | 217.9751 |
| VaR.ew.optimal | 194.3892 |
| VaR.ew.optimal.cac | 203.7821 |
| VaR.ewma | 470.9204 |
| VaR.ewma.martin | 471.0252 |
Second, the EWMA case.
Pvar2_EWMA <- t(alpha2) %*% C %*% alpha2
Psd2_EWMA <- Pvar2_EWMA^0.5
VaR.ewma.optimal <- qnorm(0.99, 0, 1) * Psd2_EWMA
VaR.ewma.optimal## [,1]
## [1,] 338.2937Pvar3_EWMA <- t(alpha3) %*% C %*% alpha3
Psd3_EWMA <- Pvar3_EWMA^0.5
VaR.ewma.optimal.cac <- qnorm(0.99, 0, 1) * Psd3_EWMA
VaR.ewma.optimal.cac## [,1]
## [1,] 347.9537Summary of results.
Results <- data.frame(VaR.method=c("VaR.hist",
"VaR.hist.quantile",
"VaR.ew",
"VaR.ew.optimal",
"VaR.ewma",
"VaR.ewma.optimal",
"VaR.ewma.optimal.cac",
"VaR.ewma.martin"),
value=c(VaR.hist,
VaR.hist.quantile,
VaR.ew,
VaR.ew.optimal,
VaR.ewma,
VaR.ewma.optimal,
VaR.ewma.optimal.cac,
VaR.ewma.martin))
kable(Results, caption = "Summary of results.", digits = 4) |>
kable_styling(latex_options = "HOLD_position")| VaR.method | value |
|---|---|
| VaR.hist | 253.3850 |
| VaR.hist.quantile | 218.3281 |
| VaR.ew | 217.9751 |
| VaR.ew.optimal | 194.3892 |
| VaR.ewma | 470.9204 |
| VaR.ewma.optimal | 338.2937 |
| VaR.ewma.optimal.cac | 347.9537 |
| VaR.ewma.martin | 471.0252 |
And third, my EWMA version.
# My EWMA approach
VaR.ewma.martin## [,1]
## [1,] 471.0252Pvar_EWMA12 <- t(alpha2) %*% VCV %*% alpha2
Psd_EWMA12 <- Pvar_EWMA12^0.5
VaR.ewma.martin.optimal <- qnorm(0.99, 0, 1) * Psd_EWMA12
VaR.ewma.martin.optimal## [,1]
## [1,] 338.4803Pvar_EWMA13 <- t(alpha3) %*% VCV %*% alpha3
Psd_EWMA13 <- Pvar_EWMA13^0.5
VaR.ewma.martin.optimal.cac <- qnorm(0.99, 0, 1) * Psd_EWMA13
VaR.ewma.martin.optimal.cac## [,1]
## [1,] 348.1087Summary of results.
Results <- data.frame(VaR.method=c("VaR.hist",
"VaR.hist.quantile",
"VaR.ew",
"VaR.ew.optimal",
"VaR.ewma",
"VaR.ewma.optimal",
"VaR.ewma.martin",
"VaR.ewma.martin.optimal",
"VaR.ewma.martin.optimal.cac"),
value=c(VaR.hist,
VaR.hist.quantile,
VaR.ew,VaR.ew.optimal,
VaR.ewma,
VaR.ewma.optimal,
VaR.ewma.martin,
VaR.ewma.martin.optimal,
VaR.ewma.martin.optimal.cac))
kable(Results, caption = "Summary of results.", digits = 4) |>
kable_styling(latex_options = "HOLD_position")| VaR.method | value |
|---|---|
| VaR.hist | 253.3850 |
| VaR.hist.quantile | 218.3281 |
| VaR.ew | 217.9751 |
| VaR.ew.optimal | 194.3892 |
| VaR.ewma | 470.9204 |
| VaR.ewma.optimal | 338.2937 |
| VaR.ewma.martin | 471.0252 |
| VaR.ewma.martin.optimal | 338.4803 |
| VaR.ewma.martin.optimal.cac | 348.1087 |
5.5 Historic approach with optimal weights.
Let’s go back to the historic approach. Now we can incorporate a comparison with the optimal weights.
# Miniumu variance optimal weights.
p2 <- (alpha2[1] * DJIA[2:501, ] / DJIA[1:500, ]) +
(alpha2[2] * AFTSE100[2:501, ] / AFTSE100[1:500, ]) +
(alpha2[3] * ACAC40[2:501, ] / ACAC40[1:500, ]) +
(alpha2[4] * ANikkei[2:501, ] / ANikkei[1:500, ])
# Optimal portfolio targeting the CAC 40 expected return.
p3 <- (alpha3[1] * DJIA[2:501, ] / DJIA[1:500, ]) +
(alpha3[2] * AFTSE100[2:501, ] / AFTSE100[1:500, ])+
(alpha3[3] * ACAC40[2:501, ] / ACAC40[1:500, ]) +
(alpha3[4] * ANikkei[2:501, ] / ANikkei[1:500, ])
# Losses in both cases.
l2 <- 10000 - p2
l3 <- 10000 - p3This is the Hull original case.
plot(l, type = "h", ylim = c(-600, 600), lwd = 2,
xlab = "500 scenarios (sorted by historical date)",
ylab = "Gains (-) and losses (+)", col = ifelse(l < 0, "blue", "red"))
abline(0, 0)
abline(v = 0)
abline(h = max(l), lty = 2)
abline(h = min(l), lty = 2)
Figure 5.14: Losses with Hull’s suboptimal alpha DJ=0.4, FTSE=0.3, CAC=0.1, Nikkei=0.2.
The one-day 99% value at risk can be estimated as the fifth-worst loss.
sort(l)[496]## [1] 253.385Now, let’s see the minimum variance optimal portfolio.
plot(l2, type = "h", ylim = c(-600, 600), lwd = 2,
xlab = "500 scenarios (sorted by historical date)",
ylab = "Gains (-) and losses (+)", col = ifelse(l2 < 0, "blue", "red"))
abline(0, 0)
abline(v = 0)
abline(h = max(l2), lty = 2)
abline(h = min(l2), lty = 2)
Figure 5.15: Losses with optimal alpha DJ=0.56, FTSE=0.038, CAC=0.022, Nikkei=0.38.
The one-day 99% value at risk can be estimated as the fifth-worst loss.
sort(l2)[496]## [1] 213.6898See how the range of losses and gains is now reduced.
And this is the optimal portfolio targeting the CAC 40 expected return.
plot(l3, type = "h", ylim = c(-600, 600), lwd = 2,
xlab = "500 scenarios (sorted by historical date)",
ylab = "Gains (-) and losses (+)", col = ifelse(l3 < 0, "blue", "red"))
abline(0, 0)
abline(v = 0)
abline(h = max(l3), lty = 2)
abline(h = min(l3), lty = 2)
Figure 5.16: Losses with CAC40 optimal alpha DJ=0.607, FTSE=-0.403, CAC=0.459, Nikkei=0.336.
The one-day 99% value at risk can be estimated as the fifth-worst loss.
sort(l3)[496]## [1] 248.8909Note that the range is reduced but allows for a higher gains. Let’s see these three together.
par(mfrow=c(1, 3), oma = c(0, 0, 2, 0))
par(pty = "s")
plot(l, type = "h", ylim = c(-600, 600), lwd = 2, xlab = "",
main = "Hull original case.", ylab = "Gains (-) and losses (+)",
col = ifelse(l < 0, "blue", "red"))
abline(0, 0)
abline(v = 0)
abline(h = max(l), lty = 2)
abline(h = min(l), lty = 2)
plot(l2, type = "h", ylim = c(-600, 600), lwd = 2,
main = "Minimum variance.",
xlab = "500 scenarios sorted by historical date.", ylab = "",
col = ifelse(l2 < 0, "blue", "red"))
abline(0, 0)
abline(v = 0)
abline(h = max(l2), lty = 2)
abline(h = min(l2), lty = 2)
plot(l3, type = "h", ylim = c(-600, 600), lwd = 2, xlab = "",
main = "Optimal targeting CAC 40.", ylab = "",
col = ifelse(l3 < 0, "blue", "red"))
abline(0, 0)
abline(v = 0)
abline(h = max(l3), lty = 2)
abline(h = min(l3), lty = 2)
Figure 5.17: Hull’s suboptimal alpha (left); optimal alpha (center); optimal alpha CAC40 (right).
A quantile comparison.
lq <- quantile(l, c(0.25, 0.5, 0.75))
lq2 <- quantile(l2, c(0.25, 0.5, 0.75))
lq3 <- quantile(l3, c(0.25, 0.5, 0.75))
ll <- data.frame(lq, lq2, lq3)
colours <- c("red", "orange", "blue")barplot(as.matrix(t(ll)), beside = TRUE, col = colours, ylim = c(-60, 60),
xlab = "Quantiles", ylab = "Gains (-) and Losses (+)")
legend("topleft", c("Hull original","Optimal minvar",
"Replicating CAC40 ER"), cex = 1.3, bty = "n",
fill = colours)
Figure 5.18: Quantile comparison of losses distribution.
This document took 4.2 minutes to compile in Rmarkdown, R version 4.2.1 (2022-06-23 ucrt).